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Abstract 


The development of machine vision based pilot aids to help reduce night ap- 
proach and landing accidents is explored in this report. The techniques developed 
in this report are motivated by the desire to use the available information sources 
for navigation such as the airport lighting layout, attitude sensors and Global Posi- 
tioning System to derive more precise aircraft position and orientation information. 
The fact that airport lighting geometry in known and that images of airport light- 
ing can be acquired by the camera, has lead to the synthesis of machine vision 
based algorithms for runway relative aircraft position and orientation estimation. 

The main contribution of this research is the synthesis of seven navigation 
algorithms based on two broad families of solutions. The first family of solution 
methods consist of techniques that reconstruct the airport lighting layout from the 
camera image and then estimate the aircraft position components by comparing 
the reconstructed lighting layout geometry with the known model of the airport 
lighting layout geometry. The second family of methods comprises of techniques 
that synthesize the image of the airport lighting layout using a camera model 
and estimate the aircraft position and orientation by comparing this image with 
the actual image of the airport lighting acquired by the camera. Algorithms I 
through IV belong to the first family of solutions while Algorithms V through 
VII belong to the second family of solutions. Algorithms I and II are parameter 
optimization methods, Algorithms III and IV are feature correspondence methods 
and Algorithms V through VII are Kalman filter centered algorithms. In order to 
take advantage of the aircraft dynamics and the multiple images available along 
the glide path, the position estimates provided by Algorithms I through IV are 


used for driving a six-state Kalman filter for providing improved estimates of the 
aircraft position and inertial velocity components. Algorithms V through VII are 
Kalman filter centered algorithms and are designed to implicitly utilize the aircraft 
dynamics and the multiple images available along the glide path. Additionally, 
Algorithm VI integrates the position information derived from a Global Positioning 
System receiver. 

Results of computer simulation are presented to demonstrate the performance 
of all the seven algorithms developed in this report. It is shown that all the algo- 
rithms meet some or all of the Federal Aviation Administration specified navigation 
accuracy requirements for various landing categories. These results show that it 
is feasible to design an accurate machine vision based night landing aid with the 
currently available technology. 
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Chapter 1 


The Need for Pilot Aids 


Landing is one of the most demanding flight regimes in fixed-wing aircraft 
operations. This fact is borne out by the statistic that the landing phase of flight 
alone accounts for 29% of all the aviation accidents. Approach and landing acci- 
dents together account for 41% of all aircraft accidents [5]. Research shows that 
night approach accident rates are about eight times that of the day rate [9]. This 
is perhaps attributable to difficulties associated with reduced lighting during the 
nighttime hours. Clearly, out-of-the-window references, navigation aids, and air 
traffic awareness are significantly impacted during these low-light conditions. Fur- 
thermore, the human body is primarily adapted for daytime activity. Night flying 
places the pilot’s eyes, which provide the primary sensory information needed for 
flight, in an environment for which they are only partially suited. Limitations of 
the human visual system along with aircraft motion are responsible for numer- 
ous static and dynamic illusions which can have dangerous consequences on night 
landing [40]. Thus, in addition to the usual landing hazards such as winds aloft, 
and complex approach procedures employed at airports near population centers, 
night landing can further add to the pilot work load. 

Landing aids such as the Instrument Landing System (ILS) and the Microwave 
Landing System (MLS) can be used to ameliorate the night landing difficulties. 
Due to their high cost, these systems are likely to be available only at a few 
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major airports. Given the operational advantages of all weather landing at any 
airport, large commercial carriers are likely to equip their airplanes with such 
systems. Currently, ILS systems are routinely used by these air carriers to land 
their airplanes. Emerging Global Positioning System (GPS) technologies hold the 
promise for low-cost, high precision landing guidance. GPS-based landing systems 
are likely to find widespread applications in the aeronautical operations. 

Smaller air carriers and general aviation aircraft which are not equipped with 
INS can only navigate along the Victor Airways or Jet Routes to the destina- 
tion airport using very-high-frequency omnirange equipment (VOR) and distance 
measuring equipment (DME) [49]. Once the destination airport is visible, runway 
lighting is used for obtaining alignment guidance. Visual approach slope indicator 
(VASI) or precision approach path indicator (PAPI) lights are used for obtaining 
glide slope information. The objective of the research given in this report is to 
explore the development of a pilot aid that can help reduce night approach and 
landing accidents. The research focus is on developing an onboard instrument that 
complements existing cockpit instrumentation. 

The techniques developed during the course of this research are motivated by 
the desire to use the emerging machine vision techniques along with the existing 
infrastructure to derive more precise aircraft state information. Decreasing costs 
of machine vision systems and components places this technology in an attractive 
position. Even if a highly sophisticated landing system were to become available, 
runway lighting will continue to be in use. Thus, the machine vision based system 
will be the ultimate back-up landing system. As and when GPS becomes cheaper 
and more accurate, the machine vision system can be used to further add value 
to it. Finally, even though the focus of this report is general aviation application, 
there is no reason why the algorithms and methods proposed here cannot be used 
in commercial and military aircraft. 

In order to further motivate the development of landing aids, factors that 
make night landing hazardous are next examined. 
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1.1 Pilot’s Health Condition 

The pilot is required to be in good health in order to cope with all the situa- 
tions encountered during night flying. The following factors are indicated in [40] 
as symptoms of changing health. Sleeping problems, chronic fatigue, gastric dis- 
turbances, shortness of breath, appetite changes, reduced eye-hand coordination 
or muscle tremor, high blood pressure, and body weight change of more than ten 
percent when not dieting. Of these, the important ones are sleep disorders and 
fatigue. 

Pilots, like other human beings, experience regular sleep and wakefulness cy- 
cles in consonance with the day-night cycles. This is known as circadian rhythm. 
This rhythm resets the biological activities once every cycle. Pilots are required 
to stay awake during night flight which conflicts directly with their need to sleep 
during the nighttime hours. Lack of sleep causes sleep disorders and fatigue. Sleep 
disorders are also caused if one’s sleep hours are shifted to a new time period dur- 
ing the day. For example, transmeridian flights require synchronization of body 
rhythms to new time zones. Usually, this adjustment is accompanied by loss of ap- 
petite and tiredness. Other side effects of sleep deprivation are short-term memory 
loss, reduced attention span, reduced judgement capability, increased irritability 
and anxiety, and increased reaction time. 

Fatigue can be defined as a general loss of well-being caused by physiological 
and psychological factors such as inadequate rest or sleep, intense mental activity, 
limited visibility, seat discomfort, airplane vibration and noise, and excessive radio 
communications. Pilot response to fatigue is very similar to that caused by sleep 
deprivation. 

Sleep deprivation, fatigue and a number of other factors that impair pilot 
performance and the physiological and psychological responses to these causal 
factors are described at length in [40]. Steps needed for preventing and overcoming 
night pilot health related difficulties are also listed in Reference [40]. 
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1.2 Flight Situational Awareness 

In order to maintain flight safety, a correct assessment of aircraft attitude is 
needed at all times. Although all aircraft include cockpit instrumentation needed 
for safe flight operations, to a large extent pilots base their sense of orientation 
on visual, vestibular and somatosensory systems. Often these reflexes give a false 
sense of attitude. Therefore, a trained pilot consciously suppresses the unwanted 
vestibular and somatosensory reflexes, and uses only the information that is visu- 
ally derived [24]. However, a number of potentially dangerous situations may be 
attributed to the information provided by vestibular and somatosensory systems. 

The vestibular system consists of the semicircular canals and the otolith organs 
of the middle ear. The semicircular canals and otolith organs provide information 
about angular and linear accelerations, respectively. In addition, the otolith or- 
gans also sense the direction of the gravity vector. The information provided by 
the vestibular system is needed for stabilization of the eyes during head or body 
motions, which would otherwise result in blurred vision. In the absence of vision, 
accurate motion and orientation perception when on ground is also provided by 
the vestibular system. Although, the vestibular system is ideally suited for the 
ground environment, it is only partially suited for the flight environment. Under 
certain flight conditions it can generate false perceptions. 

The somatosensory system responds to pressure and stretch. It consists of 
somatosensory sensors that are distributed in several body structures, including, 
skin, joints and muscles. This system is responsible for the so-called “seat— of— the— 
pants sense referred to by pilots [24]. Like the vestibular system, the somatosen- 
sory system can also generate false perceptions under certain conditions. 

In addition to the vision, vestibular and somatosensory systems, pilots learn 
to use the auditory system to get a sense of airspeed and attitude based on the 
wind noise in the cockpit [24]. 

Compared to the vestibular and somatosensory systems, the visual system 
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provides more accurate orientation information. In situations such as nighttime 
flight operations, the visual information is considerably degraded, forcing the pilots 
to depend on less accurate vestibular and somatosensory systems. In the next 
section, the impact of reduced lighting flight operations on the visual system is 
examined. 


1.3 Vision at Night 

A pilot’s vision provides the primary sensory information required for flight. 
Hence, it is important to examine how the human visual system is impacted dur- 
ing the twilight and nighttime hours. The visual information is combined with 
other sensory information, memory and domain knowledge via complex mental 
processes to result in visual perception or understanding of the scene. For cor- 
rect interpretation of the flight situation, both static and dynamic visual cues are 
needed. 

The cockpit layout around the night aviator plays an important role in pro- 
viding a frame of reference for the pilot. It is with this reference that the pilot 
perceives himself or herself to be a fixed component of the aircraft. Static structure 
of the cockpit aids the pilot in making appropriate control inputs by providing a 
stable visual reference for judging changes relative to the external environment. 
Static cockpit structure is so important that excessive head motions have been 
known to result in a sense of uncertainty about the aircraft attitude. 

Static structures external to the cockpit such as the aircraft nose restrict ex- 
ternal visibility. To overcome this difficulty to some extent, a design eyepoint is 
specified for the cockpit to permit optimal internal and external visibility. Pilots 
are required to be positioned correctly with respect to the design eyepoint. How- 
ever, over time pilots may have a tendency to slump down into the seat, thus 
lowering the eye position by a couple of inches, thereby causing a significant devi- 
ation from the design eyepoint. This is very significant during night landing since, 
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deviations from the design eyepoint can result in diminished visual range. This 
could cause the runway lights during the final approach to appear later than they 
would have if the visual range were greater. 

Spatial reference is also established by the ground plane which provides the 
horizon. Objects of known size on the ground provide scale and distance infor- 
mation. The motion of the objects in the visual field provides information about 
groundspeed. During nighttime the horizon and the objects are difficult to see. In 
some situations this can lead to a complete loss of spatial orientation. Such disori- 
entation causes symptoms of fright, airsickness and dizziness. The recommended 
procedure in such situations is to switch the pilot’s attention to the cockpit instru- 
ments. 

Due to the greatly reduced visual information during nighttime flight opera- 
tions, pilots are unable to compensate for perceptual disturbances. A major cause 
of perceptual disturbances is head motion. During and after rolling and pitching 
head motions, pilots have reported a feeling that the flight situation may be less 
safe and secure. This is probably due to conflicting information from vestibular 
sense organ and the visual system. Due to this reason, the head should be kept as 
motionless as possible. However, pilots do have to continually scan the external 
environment and cockpit instruments. Since body motions are deliberately carried 
out, any apparent motion of cockpit structures, such as window frames, relative 
to the external environment are attributed to the body motions. All other mo- 
tions are inferred to be due to aircraft motion. These two types of motion are not 
easily distinguishable by the night pilot because the visual cues needed for correct 
interpretation are either lacking or are considerably degraded during nighttime 
hours. 

The combination of reduced lighting, perceptual disturbances and the motion 
of the outside scene perceived by the pilot give rise to a number of potentially 
dangerous visual illusions. A few commonly encountered illusions are discussed 
next. Reference [40] discusses these in further detail. 
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1.4 Visual Illusions 

A visual illusion is a false perception of reality. Often, false perceptions are 
a consequence of logical interpretation by the observer. Visual illusions can occur 
when there is differential motion between the outside scene and the aircraft that 
is perceived within the pilot’s field-of-view. They also occur in situations when 
the outside scene moves across the pilot’s field— of— view during relatively stable 
visual fixation. The commonly known visual illusions that a night pilot is faced 
with during descent and final approach for landing are described in the following 
sections. These descriptions are primarily based on [40]. 

1.4.1 Runway Length/Width Illusion 

During the final approach to landing, pilots gauge the aircraft position with 
respect to the runway and the glide slope by how long and wide the runway appears 
from their viewing position. During the night, objects of known size and shape 
on the runway surface are not clearly visible. Therefore, the length/width illusion 
may arise because of what is observed differs from the pilot’s expectation. If the 
runway width appears to be larger, the pilot perceives the aircraft to be below 
the normal glide path. A narrower runway on the other hand gives the illusion of 
being high. The latter can cause the pilot to increase the rate of descent. Since the 
aircraft is close to the ground, by the time the pilot realizes that the aircraft will 
land short, there may not be enough lift margin to arrest the rate of descent [24]. 

1.4.2 Foreshortening Illusion 

Foreshortening illusion pertains to when the true shape of the objects such 
as terrain features appear to be more elliptical or shortened when viewed from a 
distance along the glide slope. 
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1.4.3 Sloped Runway Illusion 

If the ground surface is not level, as in the case of sloped runways, the visual 
cues effect pilot’s judgement of the aircraft altitude and glide slope. Usually, 
runways are at the same level as the surrounding terrain. Therefore, the visual 
information from the terrain reinforces the runway perception. In situations where 
the runway actually slopes with respect to the ground while pilots expect the 
runway to be level with the ground, it has been observed that steeper approaches 
are flown to upsloped runways and shallower approaches to downsloped runways. 
The illusionary condition causes the pilot to land short of the touchdown point 

on runways with upslope, and to overshoot the touchdown point on runways with 
downslope. 

1.4.4 False Horizon Illusion 

False horizon illusion mainly relates to misinterpretation of the location of the 
horizon within the field-of-view. One form of this illusion occurs when lights on 
the ground appear to merge with stars. This results in pilots placing the aircraft in 
unusual attitudes in an attempt to keep some ground lights above, having perceived 
them as stars. Another form of this illusion occurs when several lights are seen 
beyond the runway at a higher elevation. These lights give the impression of a 
horizon, prompting the pilots to fly below the glide slope. 

1.4.5 Vertical Position Illusion 

Well lighted objects or terrain features that are farther away from the pilot 
appear higher on the horizon. This may give the impression that the aircraft is 
higher on the glide slope than it actually is. This can result in a descent rate in- 
crease reaction. Vertical position illusion when combined with false horizon illusion 
leads to other illusions. One of these occurs when the pilot observes a light located 
on the ground a small distance ahead and to the side. The pilot may have to look 
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at a downward angle to observe the light relative to the wings. This gives the 
impression that the line-of-sight to the ground is more level which may cause the 
pilot to assume a nose up attitude. At low altitudes, as the pilot looks downward 
at the ground light when the horizon is invisible, a small bank angle may develop. 
In this case, the pilot is unable to perceive the development of this dangerous bank 
angle. 

1.4.6 Illusions Caused by Fog and Rain 

As pilots descend to the runway, presence of fog causes the runway lights to 
appear less bright, causing a misperception of the actual distance from the runway. 
Pilots are led to believe that the aircraft is farther away than it truly is. Refraction 
caused by heavy rain on the windshield results in ground lights to appear from an 
apparent location. This may give rise to errors in perceived altitude. Rain can also 
cause the runway to appear larger in size when compared to clear air conditions and 
can cause the horizon to appear closer. Heavy rain can often cause the complete 
disappearance of the horizon. When an approach is made through fog or haze, 
vertical visibility is better than forward visibility. This causes the ground lights 
farther ahead to appear less bright leading to the illusion that the aircraft is pitched 
up. 

1.4.7 After-Image Illusion 

A visual after-image remains when an observer views a bright light at night. 
For example, a camera flash bulb leaves a visual after-image subsequent to going 
off. This after-image results in the illusion that the environment is more static, 
hence, attitude changes are not perceived during this period. This illusion is 
encountered specially during approach to a runway, since, high intensity strobe 
lights placed along the runway approach centerline are used for indicating approach 
direction. Once a certain altitude is reached, pilots often request the control tower 
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to turnoff these lights. 

1.4.8 Ganzfeld Depth Loss Illusion 

Ganzfeld is a German word for a featureless visual scene. This illusion results 
in loss of depth perception when flying over snow fields, bodies of water or any other 
featureless terrain. Some features are required on the terrain so that location of 
one feature is judged to be different than the location of another feature. Without 
this prerequisite, depth discrimination is impossible. At night, unilluminated areas 
of the terrain with vastly different features appear continuous. For example, bodies 
of water smoothly merge with land in the visual scene. 

1.4.9 Foreground Occlusion Illusion 

This illusion is most often experienced when the ground lights are cutoff by a 
cloud. In a moonlit night, pilots can detect the cloud by its reflection. However, 
in a dark night such discrimination is not easy. A more dangerous version occurs 
during descent at night over mountainous terrain. During a portion of the descent, 
the lights on the runway are visible to the pilot and the foreground occlusion such 
as a hilltop lies invisible. At some point along the descent, the lights are cutoff 
by the hilltop. When such a situation is detected, pilot should climb immediately 
or else a collision with the terrain would occur. It is easy to see how this illusion 
could cause confusion in judging whether a hilltop or a cloud is the cause for the 
foreground occlusion. Detailed terrain knowledge is one of the useful sources of 
information for correct interpretation. 

1.4.10 Up— Sloped Lighted City Illusion 

This illusion is experienced when terrain stays level for some distance and then 
rises to give the impression of two intersecting planes in the pilot’s field-of-view. 
Often there are parallel roads with street lights in a city situated on the upward 
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sloping terrain. Long rows of street lights appear to converge at a distance giving 
the impression of a horizon. The runway lights, situated on the level terrain in the 
foreground, may also appear to converge at a different vanishing point. The pilot 
can be tricked into believing that the broad upward sloping terrain is level and 
that the runway is sloped down. This may cause the pilot to increase the descent 
rate. 

1.4.11 Autokinetic Illusion 

When an observer stares steadily at a single motionless source of light at 
night, autokinetic illusion gives the appearance that the source of light is moving 
around in random directions at varying speeds. Due to this illusion, an isolated 
motionless ground light may appear to be moving on the ground. One possible 
erroneous interpretation is that another aircraft is in the vicinity. Autokinetic 
illusion can also cause a visible star to be misperceived as a moving vehicle on the 
ground, giving the impression of low pitch attitude to the pilot. 

1.4.12 Black Hole Approach Illusion 

Black hole approach illusion arises during night approaches where no ground 
details are visible short of the runway. Four different types of black hole approach 
situations have been described in [40]. The main factor that causes this illusion is 
that pilots derive vertical guidance information in the angle between the line-of- 
sights to the farthest and the nearest light. If an aircraft is flown at a constant 
altitude, the angle is expected to increase as the aircraft nears the runway. Simi- 
larly, the angle should decrease as the aircraft descends. In cases where the pilot 
is unable to perceive visual angle change, a more rapid descent is flown. Problem 
occurs in situations where the aircraft descends into the terrain much before the 
runway. 

This concludes the discussion of visual illusions. A few vestibular and so- 
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matosensory illusions are described next. 


1.5 Vestibular and Somatosensory Illusions 

Vestibular and somatosensory illusions are caused by the linear and angular 
accelerations perceived by the pilot. 

Somatogravic illusion is a sensation of change of attitude when the otolith 
organs are subjected to linear acceleration. This illusion occurs in level flight 
giving the pilot a false cue of being in a nose high attitude during acceleration. 
The opposite illusion of nose down attitude occurs during deceleration on final 
approach. A pilot may create a low altitude stall in the process of correcting for 
this illusion [24]. A variant of this illusion is the inversion illusion in which the G 
forces acting on the otolith organs give the sensation of being upside down, when 
the pilot is being subject to negative G forces [24]. 

During a coordinated turn, the “seat-of-the-pants” sense is misleading be- 
cause the resultant of the gravitational and centrifugal forces is directed towards 
the floor of the aircraft, which the pilot falsely perceives as the direction of the 
vertical [24]. 

Coriolis illusion occurs during prolonged turns in one plane. The sensation 
of turn perceived by the semicircular canals in the inner ear at the beginning of 
the turn subside during the prolonged turn. A sudden head movement causes 
the canals to sense angular acceleration which gives the impression of rotation 
in another plane. Attempts to correct for this illusion can place the aircraft in 
dangerous attitudes [24]. The coriolis illusion is specially hazardous during curved 
approach because of the aircraft’s proximity to ground. Furthermore, it can cause 
disorientation and can produce intense symptoms of nausea [40]. 

“Leans” is a common illusion caused by rapid roll maneuvers to correct for roll 
angle developed by subtle bank. For example, if a subtle bank angle develops to 
the left such that the vestibular system is unable to detect it, the pilot eventually 
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detects the roll on the attitude indicator and corrects it by a rapid roll to the right. 
The pilot gets the false sense of only having rolled to the right. 

Giant hand illusion gives the impression that an external force is pushing on 
the aircraft or holding it at a certain attitude. This is caused by vestibular and 
somatosensory inputs that interfere with pilot’s conscious control of the aircraft. If 
the disorientation is about the pitch axis during aircraft acceleration, the aircraft 
appears to resist pilot efforts to pull the nose up because the natural reflex is to 
push the nose down [24]. This illusion also occurs when the disorientation is about 
the roll axis as in the “Leans” illusion. In these cases, the aircraft seems to resist 
roll efforts by the pilot. 

In addition to the illusions described here, a number of vestibular and so- 
matosensory illusions can occur in high performance aircraft during maneuvers 
such as graveyard spin, graveyard spiral and rapid aileron rolls. These are de- 
scribed in further detail in Reference [24]. 

With the background of sensory illusions that the night pilots often experience, 
the potential use of machine vision systems in ameliorating the impact of these 
illusions is examined next. 


1.6 Machine Vision Systems As Pilot Aids 

Based on the preceding discussion of the human perceptual system and how 
prone it is to visual, vestibular and somatosensory illusions, this report attempts to 
answer the question: Can a machine vision system augment the pilot’s perception 
sufficiently to avoid these illusions? 

Before an attempt is made to answer this question, it is necessary to establish 
the underlying causes for the various illusions described earlier. Closer analysis 
reveals that they can be classified into three distinct groups based on the underly- 
ing causal factors. Those that occur because of imprecise knowledge of geometry, 
those due to conflicting information from the vestibular, somatosensory and vi- 
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sual systems, and those arising from the limitations of the human eye. Runway 
length/width illusion, foreshortening illusion, sloped runway illusion, false horizon 
illusion, vertical position illusion, Ganzfeld depth loss illusion, foreground occlu- 
sion illusion, up-sloped lighted city illusion, and black hole approach illusion, all 
have their root in imprecise knowledge of the terrain, runway and lighting geom- 
etry. Autokinetic illusion and the various vestibular and somatosensory illusions 
have their roots in conflicting information received from the vestibular, somatosen- 
sory and vision systems. The third group of illusions, which includes after-image 
illusion and fog and rain caused illusions, has its basis in the physical limitations 
of the human eye. After-image illusion is caused due to saturation of the visual 
receptors in the eye. Fog and rain cause the runway lights to appear diffused. In 
this situation, the eye has no mechanism for enhancing the appearance of these 
lights. 

In order to examine how a machine vision system could be functionally supe- 
rior to its human counterpart, consider the analogy between the human perceptual 
system and the machine vision system shown in Figures 1.1 and 1.2. 




EYE 



VESTIBULAR ORGAN 


Figure 1.1: Human perceptual system. 


The human perceptual system is mainly driven by three sources for the land- 
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Figure 1.2: Proposed machine vision system. 

ing task: runway lighting as seen by the eye, motion sensed by the vestibular and 
somatosensory systems, and runway knowledge learnt and stored in the memory. 
An analogous machine vision system could also be driven by equivalent sources: 
runway lighting seen by the camera, motion sensed by accelerometers and gyro- 
scopes, and knowledge of runway geometry encoded in the computer memory. In 
the human perception system the brain integrates the input information. In the 
machine vision system this function can be accomplished by computer-resident 
algorithms. Thus, according to this analogy, the camera can be considered equiva- 
lent to the human eye, accelerometers and gyros equivalent to the vestibular organ, 
and the geometry information available from digital memory can be considered to 
be equivalent to the domain knowledge in the human brain. 

The exact geometric information encoded in the digital memory of the com- 
puter is precise when compared with the approximate runway geometry knowledge 
stored in the human brain. For the machine vision system, this fact offers the im- 
munity to visual illusions caused by imprecise knowledge of the runway geometry. 

Accelerometers and gyros are precision instruments which far exceed the ca- 
pabilities of the human vestibular and somatosensory systems. In addition, these 
sensors provide true motion of the aircraft. In situations where the pilot moves 
relative to the airplane, the vestibular and somatosensory systems sense a combi- 





16 


nation of the aircraft motion and the pilot’s motion. Thus, it is difficult for the 
pilot to differentiate between self motion and aircraft motion. 

Unlike the human eye, the camera sensor elements can be designed to have 
optimal sensitivity to runway lights. Additionally, optical filters can be used for 
reducing or eliminating certain frequencies from the visible spectrum. They can 
also be designed to avoid saturation of the photosensitive elements. Thus, after- 
image illusion can be effectively eliminated in a machine vision system. 

So far, information sources which have the potential of providing superior 
quality information to the machine vision system have been discussed. However, 
the critical component of a machine vision system is the algorithm for estimating 
runway relative position and attitude of the aircraft. The point of view adopted 
in this research is that two categories of algorithms based on sound physical and 
mathematical principles are needed for algorithm development. Firstly, methods 
for conditioning the image output from the camera are required. Secondly, meth- 
ods for integrating information from the image, motion sensors, and the stored 
runway geometry, for runway relative position and attitude determination need 
to be developed. Both categories are topics in the Computer Vision or Machine 
Vision literature. Assuming that such algorithms can be designed, the next issue 
relates to what is available in the literature and what has been accomplished so 
far. The following chapter provides a brief summary of Computer Vision and the 
literature relevant to the design of such algorithms. 

The preliminary discussion in this section provides a glimpse at the possibil- 
ities offered by a machine vision system. Although current generation imaging 
sensor technology is adequate for the design of a machine vision system, future 
improvements will only enhance the capability of such a system. 
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1.7 Summary 

The complexity of the landing task and the hazards of night operations were 
discussed in this chapter. Pilot health issues, flight situational awareness, vision 
at night, visual illusions, vestibular and somatosensory illusions were discussed at 
length. A study of these issues indicated that a low cost, machine vision based 
position and orientation instrument was required for general aviation. Since human 
pilots are able to fly the aircraft along the descent path purely by visual reference, 
it was argued by drawing an analogy between the human perceptual system and the 
machine vision system that a machine vision system could be designed for deriving 
runway relative position information during approach and landing without being 
subject to optical illusions. Reasons were given for expecting higher reliability of 
the machine vision system as opposed to the human perceptual system, specially 
in cases where precise knowledge of runway geometry is required. It was pointed 
out that the algorithms are the key component of the machine vision system. 



Chapter 2 


The Machine Vision Technology 


Machine or Computer Vision technology deals with algorithms and methods 
for two dimensional digital signal processing, pattern classification, image segmen- 
tation, geometric modeling, and relational structures. Text books in this area 
[6, 27, 47, 53, 78] cover many of the topics of machine vision. Many of these text 
books have an Artificial Intelligence flavor, focusing on heuristics of machine vi- 
sion technology, Reference [6] being an example of this approach. A few examine 
the issues from a signal processing point of view, while others are motivated by 
statistical decision theory. Representative examples of these two approaches are 
References [53] and [27]. These references cover most of the tools and techniques 
used in machine vision system development and research. 

The range of topics addressed in machine vision technology can be organized 
as a sequence of representations from early or low-level vision to cognitive inter- 
pretation [6]. Starting with an image generated by the camera, early or low-level 
vision algorithms are used for image conditioning such as filtering, edge detection, 
and optical flow computation. The output of this process is usually encoded into 
an image format, often called intrinsic or generalized image [6]. Higher level al- 
gorithms use these intrinsic images as inputs and gather regions within the image 
that are likely to be associated with objects being viewed. For example, neigh- 
boring pixels in the image which have the same color can be grouped together to 
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represent an object. Higher-level algorithms also assign geometric properties such 
as shape, area, eccentricity, compactness, and boundary length to regions within 
the image. This representation is natural if a database containing all the shapes 
making up the scene are available. Shape properties can be used for matching the 
image with models in the database, permitting the derivation of the location and 
orientation of the observer. Finally, higher level algorithms may use rules of logic 
to infer about what is seen. Clearly, this function is very much dependent on the 
domain. With this brief introduction to machine vision, modern solid state imag- 
ing systems are examined next. Some of the low-level and high-level functions are 
examined. 


2.1 Modern Solid State Imaging Systems 

Electronic imaging technology has changed considerably since the introduc- 
tion of photoemissive storage tubes which use incident light to emmit electrons 
in a pattern corresponding to the brightness of the scene. The Iconoscope was 
the first practical device of this type. This was soon replaced by Image Orthicon. 
The low signal-to-noise ratio of these devices led to the development of photo- 
conductive devices. Photoconductive devices are based on principle of change in 
electrical resistance of a photoconductor when exposed to light. Vidicon, Plumbi- 
con and later Saticon are devices of this type. More recently, solid-state devices 
called charge-coupled devices (CCDs) have found an increased use in the consumer 
electronics. These devices provide good signal-to-noise ratio along with the ad- 
vantages of small size, low power consumption and low cost. A modern solid state 
CCD camera unit is shown in Figure 2.1. 

Cameras convert electro-magnetic radiation received within a certain field-of- 
view into electrical signals encoded to form a two dimensional array. This general 
definition is applicable to visible-range and infra-red camera systems. Thus, one 
way to classify an imaging system is by its operating range within the electro- 
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Figure 2.1: A Modern Solid State Camera. 

magnetic spectrum. This report will be mainly concerned with image sequences 
generated by Television (TV) cameras, although some of the algorithms are appli- 
cable to the infra-red systems as well. 

Major components of a solid state TV camera such as the one shown in Figure 
2.1 are the lens, iris, shutter, photosensitive sensor array and camera electronics. 
In addition to these, color cameras employ beam splitters. A compound lens is 
used for adjusting the focal length for projecting the image of the viewed scene on 
to the photosensitive sensor array. Iris controls the amount of light that is allowed 
to reach the photosensitive sensor array. Photosensitive sensor array is the sensing 
element that converts images into electrical signals. Imaging systems are often clas- 
sified by the type of photosensitive sensor used, charge coupled devices and charge 
induced devices (CID) being two examples. The camera electronics provides timing 
signals for shuttering and downloading the signal from the photosensitive sensor 
array, noise removal, signal conditioning, pre— amplification, amplification, image 
encoding and several other functions that are needed for generating acceptable 


images. 
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A CCD sensor array consists of rows and columns of photosensitive elements 
arranged in a rectangular array on a silicon substrate. Pixel size is the term used for 
describing the size of an individual photosensitive element. These elements collect 
and store electrical charge in proportion to the intensity of the light incident on 
their surface. The charges are electronically transferred to the device output to 
form the output video signal. 

The resolution of a CCD camera depends on the number of photosensitive 
elements. The number of rows determine the vertical resolution and the number of 
columns determine the horizontal resolution of the camera. Typically, a National 
Television Systems Committee (NTSC) format CCD camera is designed with 492 
rows and 510 columns [52], 

Three commonly used architectures for CCD transfer and readout are: (1) 
frame transfer, (2) interline transfer and (3) frame and interline transfer [52]. 
Frame transfer architecture uses a sensor array, a storage register array and a 
horizontal output register. The sensor array is allowed to collect charge for a com- 
plete frame. Commanded by a clock, the charges in each column are transferred 
vertically from element to element to a corresponding column of the storage array. 
This process empties the sensor array for the next frame. The charges from the 
storage array are transferred one row at a time to the output register in synchro- 
nism with clock commands. The output register generates the video signal. 

The main disadvantage of this architecture is that the photodiodes saturate if 
exposed to bright light for a long duration. To overcome this problem, a mechanical 
shutter is employed to permit the CCD sensor array to be exposed to bright light for 
a short duration. Clearly, this introduces a mechanical element into an otherwise 
all electro-optic device. This is specially of concern if the camera is to be used in 
a high vibration environment. 

In the interline architecture, every photosensitive element in the sensor array 
is connected to a neighboring storage element. The storage elements are arranged 
in columns next to sensor element columns. Once charge is collected, commanded 
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by a clock, charge is transferred horizontally from each sensor element to the 
neighboring storage element which frees up the sensor element for the next frame. 
Charge is then transferred vertically from storage element to storage element in 
each storage column. Finally, under a clock command, the charges from the storage 
array are shifted one row at a time into the horizontal output register, similar to 
the frame transfer architecture discussed earlier. The video signal is then read out 
from the output register. 

This architecture has the advantage of being resistant to blooming and smear 
because of the rapid transfer from the sensor array to the storage array. The main 
disadvantage is that placement of storage elements next to sensor elements causes 
the sensor to have lower pixel density. 

Frame interline architecture employs a row of selection gates between the 
sensor and storage elements so that excess charges are drained from the system 
before being transferred to the storage columns. This architecture is similar to 
the frame transfer architecture with the added advantage of being resistant to 
blooming and smear. 

Blooming occurs when a CCD sensor element saturates and spills charge to 
the neighboring elements. This gives the appearance of a large bright spot in the 
image. This effect may be seen in the image of the runway lighting, acquired by 
a CCD camera, shown in Figure 2.2. In order to overcome this problem, more 
expensive CCD sensor elements are designed with built in anti-blooming gates 
which remove the excess charge. 

Sensitivity is an important measure of camera performance. Camera sensitiv- 
ity is defined as the amount of light that is needed to produce a video signal of 
certain magnitude. For example, sensitivity can be characterized by the amount 
of light in units of Candle Power required to produce a gray-level of 255 in an 
8 bit system. A more sensitive camera requires less amount of light to produce 
the same output as a less sensitive camera. Camera sensitivity can be adjusted 
to a certain extent by increasing the video gain in the camera. The disadvantage 
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Figure 2.2: A CCD camera image during night landing illustrating the “blooming” 


effect. 
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in increasing sensitivity is that it decreases the signal-to-noise ratio. For CCD 
cameras, signal-to-noise ratio is directly related to the camera sensitivity. 

Dynamic range of a CCD camera is a measure of its range of useful operation. 
It is defined as the ratio of the number of electrons required for maximum charge 
to the number of electrons that accumulate at the charge site if no light is incident 
on it. This ratio is often expressed in decibels (dB). As an example, if 80,000 
electrons represent full charge and 20 electrons represent the dark current, the 
dynamic range is 72 dB. 

Integration time is defined as the duration in which charge is allowed to ac- 
cumulate in the charge sites of the CCD array. The integration time is controlled 
by electronic shuttering or by the selection of the readout pulse width. 

Gamma is another term commonly associated with TV cameras. An image 
gamma of unity means that the system correctly reproduces the gray-levels of the 
scene [52]. If gamma is greater than unity, the image appears sharper but the scene 
contrast range is reduced. Reduction of gamma makes the image appear washed 
out [52]. Since the CCD is a nearly linear device, its output signal is directly 
proportional to the scene illumination. However, the phosphors used in display 
monitors behave nonlinearly. Typically, they have lower gain for dark signals and 
higher gain for bright signals. To compensate for this, higher gain is used for dark 
signals and lower gain for bright signals in the video output to produce a linear 
response from the monitor. This intentional incorporation of nonlinear gain is 
called gamma correction. The disadvantage of gamma correction for dark signals 
is that the noise is also amplified. Gamma correction is seldom used for cameras 
used in image processing applications. 


2.2 Low-level Vision 

The digital image generated by the camera can be considered to be a two 
dimensional function f(x,y). In order to restrict the memory requirements, it is 
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customary to represent images as an integer function of integer variables. Such 
an image with function values between 0 and 255, known as gray-levels, is shown 
in Figure 2.3. This figure is a daytime image of a runway taken by a camera 



Figure 2.3: Vehicles parked on a runway. 

mounted on the nose of a rotorcraft. Considering the image as a two-dimensional 
signal permits the application of various signal processing techniques. Techniques 
such as low-pass, band-pass and high-pass filtering, histogram equalization, and 
interpolation can all be used. Mathematical tools of transform theory such as two- 
dimensional Fourier transform, sine transform, cosine transform, singular value de- 
composition, and Radon transform can all be applied to enhance the information 
content in an image. Ideas from the theory of vector spaces can also be applied if an 
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image is conceptualized as a matrix. Images can also be processed using stochas- 
tic signal processing tools such as covariance models and autoregressive models. 
Application of some of these techniques are discussed at length in Reference [53]. 

2.2.1 Filtering 

As an example of low-level vision processing, a low-pass filtered version of 
Figure 2.3 is shown in Figure 2.4. It is hard to tell the difference between this 


Figure 2.4: Low-pass filtered image. 


image and the original image given in Figure 2.3, except for a slight reduction in 
the gray-level bandwidth. However, the difference is clear when examined in the 
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histogram domain shown in Figure 2.5. The histogram summarizes the frequency 
with which a certain gray-level appears in an image. Comparison of the two 
histograms in Figure 2.5 reveals that the gray-levels in the image shown in Figure 
2.4 vary much more smoothly. 


o 



Figure 2.5: Histograms of original image and low-pass filtered image. 


2.2.2 Edge Detection 

An important early vision processing function is the edge detection. Edges 
in an image occur at locations of large gray-level changes. These changes can be 
characterized as step, ramp or a combination of step and ramp functions. Rather 
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than attempting to locate the edges directly from the gray-level image, a gradi- 
ent operation followed by the thresholding operation is usually employed for edge 
detection. Edge operators lie in the following three classes: (1) operators that 
approximate the mathematical gradient operation, (2) template matching meth- 
ods that use many templates at different orientations, and (3) operators that use 
parametric edge models for fitting local intensities [6]. Machine vision literature 
abounds with edge operators [78]. An example of edge operation on the image in 
Figure 2.3 is given in Figure 2.6. In this case a Sobel edge operator [6] was used. 



Figure 2.6: Sobel edge magnitude for Figure 2.3. 


The Sobel edge operator is defined as: 

A u = 2(/(i + l,i)-/(i - l,j)) + (/(* + l,j - 1) -/(*- l,i - 1)) 
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+ (/(* + i,i + 1) - /(* — l, j + 1» C 2 - 1 ) 

A„ = 2(/(i,; + l)-/(i,j - 1)) + (/(i + l,i + 1) -/(* + - 1)) 

+ (/(*' - l,j + 1) - /(* - lii - 1)) C 2 - 2 ) 

with magnitude: 

s(i,j) = (A* + A*)* (2.3) 

and direction: 

x(hj) = tan_1 (^) ( 2 - 4 ) 


Here f(i,j ) is the gray-level at a pixel location (i, j). The other indices refer to the 
eight neighboring pixels surrounding this pixel. Figure 2.6 shows the thresholded 
edge magnitude. The edge direction from Equation (2.4) is illustrated as an image 
in Figure 2.7. This pseudo image is created by quantizing and scaling the edge 
direction in the range of zero and 255. Black corresponds to the vertical direction 
and white corresponds to the horizontal direction. 

2.3 Higher-level Vision 

Higher-level vision algorithms address the problems related to object repre- 
sentations in a scene. They include boundary detection, segmentation, grouping, 
geometric modeling, inference techniques and ranging. 

2.3.1 Boundary Detection Methods 

Boundary detection methods usually fall into one of the following categories: 
searching near an approximate location, Hough transform, graph searching, dy- 
namic programming and contour following [6]. There is an abundance of literature 
describing these methods [2, 8, 36, 41, 62, 65, 72]. 

Searching near an approximate location involves the determination of a likely 
location of a boundary, which is then used for guiding the refinement of the bound- 
ary. One of the methods described in Reference [6] carries out local searches at 
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Figure 2.7: Graphical representation of Sobel edge direction for Figure 2.3 
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regular intervals along directions normal to the initial boundary. In order to refine 
the boundary, the edges with the highest magnitude whose orientations are nearly 
tangential to the initial boundary are approximated using a polynomial. Recur- 
sive techniques that construct a boundary by first connecting two edges with a 
straight line and then searching along the normal at the central point for an edge 
that exceeds some preselected threshold have also been reported. This technique 
is then applied to the two segments formed by three edges and so on. Thus, a 
curved boundary is found. 

Hough transform [6, 27] can be used if little information is available about 
the location of the boundary but its shape can be described as a parametric curve. 
To illustrate the method, consider a straight line in the parameterized form: p = 
x cos 0 + y sin 0 where, 0 is the angle of the normal and p is the distance from 
the origin [27]. Given a set of points (x, •,?/,), Hough transform involves setting 
up a two-dimensional accumulator array A(9,p) which is incremented each time 
the particular (0, p) location is visited. 0 is quantized and varied between 0 and 
27 T. Hence for each (xj,y;), p’s are computed using the parameterized form and 
the accumulator array is incremented by one. If many points lie on a straight line 
corresponding to a particular 0 and p, the accumulator value for this 0 and p is 
high. Thus by using a threshold, meaningful lines in the image can be determined. 
As discussed in [6] Hough transform method can also be tailored for other shapes. 

Graph searching techniques involve searching through a set of nodes linked via 
edges to determine minimum cost paths for boundary determination. Minimum 
spanning tree algorithm described in [75] is one such graph searching algorithm. 
A spanning tree is defined as a connected set with no loops that contain all the 
points in the problem. The minimum spanning tree of a set is that whose cost 
is a minimum. On a historical note, the graph search problem is closely related 
to the travelling salesman problem in Combinatorics [59]. Several cost functions 
that can be used for boundary search are described in [6]. Heuristic graph search 
techniques and methods for managing the data structure are also described in [6]. 
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The dynamic programming procedure can also be used for boundary detection. 
This procedure builds paths from multiple starting points in a discrete region using 
a performance index that describes the optimal boundary. A recent application 
of the dynamic programming procedure for boundary detection is described in [2]. 
“Energy” is used as the performance index in that work. The energy consists 
of image intensity, edge magnitude, curvature and smoothness constraint. In this 
formulation, a penalty is imposed for moving away from the initial contour position. 

The central idea behind contour following is to start at an edge and develop 
a boundary by recursively adding neighboring edge elements based on their edge 
magnitude and direction. These methods make use of several heuristics. Recent 
methods that implement this idea are described in [8] and [62]. 

The boundary detection procedure proposed in [8] attempts to modify pa- 
rameters of lower level processes such as edge contour tracking using higher level 
processes such as corner detection. The method encodes each edge element by its 
relationship to its neighbors using a chain code scheme. A window is then used 
to determine if the neighbors extend the edge in a straight line. If the neighbors 
do not extend the edge in a straight line, left and right extensions are examined. 
This process is continued till either the contour is closed or all the pixels have been 
examined. 

A three step edge detection process is described in [62]. The first step involves 
computing the gradient magnitude and direction. The direction of the gradient 
is discretized to one of the eight neighbors surrounding the pixel. A heuristic 
concept of Likelihood— of— Being— an— Edge (LBE) is introduced as the second step. 
The third step is a contour following process which attempts to propagate the 
edge in a direction normal to the gradient direction starting at pixels with the 
maximum LBE value. The boundaries detected by application of this algorithm 
to the vehicle image is shown in Figure 2.8. 
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Figure 2.8: Edge boundaries for Figure 2.3 using the three-step process [62] 
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2.3.2 Gray-Level Segmentation 

Segmentation methods serve to consolidate the information available in the 
image. Segmentation methods either work directly with the gray-level image or 
with texture properties. The central idea is to determine contiguous regions within 
the image that satisfy some homogeneity property. Methods reported in the litera- 
ture can be broadly classified into local techniques, global techniques and splitting 
and merging techniques [6]. 

Local techniques attempt to place pixels in a region based on some homo- 
geneity property of a pixel and its neighbors. An example of this technique is blob 
coloring given in [6]. The technique involves scanning the image from left to right 
and top to bottom with a special L-shaped template. The idea is to grow the 
blob by adding a neighboring pixel if its gray-level is approximately equal to the 
gray-level of the blob. 

An example of a global segmentation technique is Thresholding. The idea here 
is that if an image consists of a background and an object, or equivalently, if the 
gray-level histogram of the image is bi-modal, a single threshold can be identified 
for segmenting the image into background and object regions. A more recent 
algorithm that extends this idea by using multiple-level thresholding is described 
in [57]. The difficulty with this technique is that many regions are given the same 
label because groups of pixels in different regions of the image lie within the same 
portion of the gray-level histogram. This is to be expected because of the global 
nature of the algorithm. Clearly, this algorithm is suitable if several objects with 
similar gray-level properties are expected to be seen against a common background 
in the image. This would be the case in nighttime images of the airport lighting. 

Segmentation methods based on merging or bottom up, splitting or top down 
and a split and merge scheme are discussed in [50]. Merging involves labeling 
of adjacent regions into a larger region based on similar gray-levels while split- 
ting involves re-labeling a larger region into several smaller regions based on the 
dissimilarity of the gray-levels in the original larger region. A split and merge 
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technique uses both splitting and merging operations. Usually merging operations 
are computationally efficient when compared with splitting operations. On the 
other hand, it requires larger memory compared to the splitting scheme. The split 
and merge scheme trades off computational speed for reduced memory requirement 
when compared to a pure merging scheme. Usually the split and merge scheme 
such as in Reference [50] is initiated at an intermediate level, close to the final 
segmentation. In this algorithm, the image is examined at various resolutions. 
Thus, four neighboring regions are merged if the difference between the maximum 
and the minimum gray-levels is less than a preset threshold. Similarly, a region in 
which the difference between the maximum and the minimum gray-levels is greater 
than a preset threshold is split into four subregions. Since both split and merge 
operations are done simultaneously, regions that are split are not merged with ad- 
jacent regions. To overcome this difficulty, a grouping technique which abandons 
the tree structure is used. Finally, the remaining small regions are merged with 
the adjacent large regions. For the vehicle image shown in Figure 2.3, 214 regions 
found by this algorithm are shown in Figure 2.9. The artifacts of the segmentation 
boundaries can be seen in this segmented image. Clearly, the power of segmenta- 
tion is also illustrated by the fact that 262144 pixel regions are compressed into 
214 regions. 

Although both boundary detection and segmentation are related, it may be 
noted that the results generated by these methods are quite different as evidenced 
by Figures 2.8 and 2.9. It is easily seen that not all boundaries are closed but all 
segmented regions are contiguous. 

The foregoing techniques are also applicable to the problem of texture segmen- 
tation with the homogeneity criteria based on texture properties. Since texture 
segmentation forms a large body of work in machine vision literature, it is treated 
here separately from gray-level segmentation. 
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Figure 2.9: Segmented regions for Figure 2.3 using the split and merge scheme 

[50]. 
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2.3.3 Texture Segmentation 

Gray-level images are characterized by pixels of varying intensity. Hence, they 
can also be described by the stochastic properties of the gray-levels distribution 
across the image. The properties of this distribution are usually given in terms 
of the first, second and higher order statistics. First order statistics describes the 
pixel population in the image without regard to its spatial distribution. The second 
order statistics take the spatial distribution into account. Two approaches are used 
to characterize this spatial distribution: (1) a stochastic model-based approach and 
(2) a data-driven approach. The model-based approach assumes that the image 
can be modeled in terms of a two-dimensional random field. Several stochastic 
models are discussed in References [43, 79]. 

The data-driven or non-parametric approach is based on characterizing the 
two-dimensional intensity distribution by different types and features of second 
order statistics. The conditional probability density function f(i,j\d,0) represents 
the probability that two pixels separated by an inter-pixel distance d and orienta- 
tion 0 have intensities i and j . An estimate of the conditional probability density 
function c(i, j\d, 0) is sometimes referred to as the gray-level co-occurrance matrix 
(GLCM) or as the spatial gray-level dependence matrix (SGLDM). SGLDM has 
been most widely used measure for classification of textures [1, 19, 39, 42, 101]. 
SGLDM can be obtained by computing the two-dimensional histogram of the 
frequency of the joint occurrences of two pixels with a fixed displacement and 
orientation with respect to each other having intensities i and j respectively. A 
rotationally invariant SGLDM is computed by averaging the individual SGLDM 
for each angular direction. 

Either matrix features or scalar features can be used for texture classification. 
Many different approaches are available for texture classification using matrix fea- 
tures. Threshold selection based on the SGLDM is described in Reference [1]. In 
Reference [19], the SGLDMs of four neighbors in the quad-tree are compared with 
a threshold for merging or splitting operations. Results using this technique are 
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also presented in [103]. A technique for image segmentation by detecting clusters 
in the SGLDM, which correspond to the regions and boundaries in the image, 
is described in [39]. A maximum likelihood texture classifier using matrix and 
scalar features is examined in [101]. In Reference [77] segmentation is carried out 
by thresholding. The threshold levels are selected by projecting the off-diagonal 
elements of the SGLDM onto the diagonal and treating the resulting vector as 
a histogram. Although these methods are useful for segmentation, their storage 
requirements are high due to the use of matrix features. For example, 256 x 256 lo- 
cations are needed to store a matrix feature for an image containing 256 gray-levels. 
These methods are also computationally intensive. The storage requirements and 
computational speed are the motivating factor for considering scalar features for 
image segmentation. However, it should be noted that many of the scalar fea- 
tures derived from the matrix features may not contain all the important texture 
information contained in the matrix features [21]. 

Several scalar features are derivable from the matrix features. For example, 
14 scalar texture features based on the SGLDM are presented in [42]. For each 
of the scalar features, their means and variance computed by using the SGLDMs 
corresponding to the four directions, may be used for texture classification. Some 
scalar features derived from SGLDM, Fourier power spectrum, Gray-level differ- 
ence statistics and Gray-level run length statistics are described in [21, 102], Scalar 
texture features derived from the SGLDM may also be computed from sum and 
difference histograms [99]. Compared to computing the full SGLDM, sum and 
difference histograms are computationally fast and require significantly reduced 
storage. Except for the two scalar features energy and entropy, all the other scalar 
features can be obtained by using the sum and difference histograms. Methods 
such as [19] and [101] can be used for classification using scalar features. Addi- 
tional methods such as piecewise linear discriminant function method, min-max 
decision rule method [42] and Fisher linear discriminant technique [102] can also 
be used for classification using scalar features. 
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Some of the scalar features relate to specific characteristics in the image such 
as homogeneity, contrast and organized structure. Other features characterize the 
complexity. Even though each scalar feature contains textural information, it is 
hard to identify which specific textural characteristic is represented by which fea- 
ture. In Reference [93], the classification characteristics of scalar features derived 
from SGLDM are examined. It is shown that scalar features used in combination 
result in superior image segmentation when compared with a single scalar feature. 
This completes the discussion of the various tools and techniques available for 
image segmentation. 

2.3.4 Clustering Methods 

The need for clustering occurs naturally in many systems. For example, vision 
based range computations [70, 94] often result in a sparse set. A vision based 
ranging method described in Reference [94] is able to compute ranges at discrete 
locations shown as white squares in Figure 2.10. Scene understanding, navigation 
and display functions require these discrete set of ranges to be grouped into sets 
which correspond to objects in the real world. Clustering techniques can be used 
for grouping the discrete range points, varying from a few hundreds to several 
thousands, into a small number of objects in the scene. 

Clustering [4] has been used for a long time in disciplines such as biology, ge- 
ology and psychiatry. In computer vision, clustering methods have been used for 
classification of multispectral data and image segmentation using attributes like 
gray-level, color, texture, gradient, and range. The main idea behind clustering or 
grouping is similar to segmentation in the sense that both the techniques attempt 
to partition a given set into subsets based on discriminants. In computer vision, 
clustering has been associated with statistical pattern recognition using discrete 
samples as in Reference [27] while segmentation has been associated with parti- 
tioning the image into homogeneous regions as in Reference [6]. In this section, 
clustering is described as a problem of partitioning discrete data with the range 
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Figure 2.10: Range locations in the image from Reference [94]. 
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map given in Figure 2.10 as an example. 

Clustering techniques can be broadly classified into supervised (model-based) 
and unsupervised (data-driven) methods. Supervised methods require labeled 
training samples. For example, if a mixture is known to be composed of samples 
from two Gaussian distributions and the problem is to separate the two types, 
known samples from each Gaussian distribution can be used for estimating the 
mean and the variance of the Gaussians. Thus, a threshold or decision bound- 
ary can be found for classifying an unknown sample into one of the two types. 
Euclidean distance and Mahalanobis distance [27] measures can also be used for 
classifying the unknown sample to the type represented by the closest mean [27]. 
In the case of unsupervised clustering, the structure is directly obtained from the 
data. However in order to design a reasonable classifier, assumptions are invari- 
ably needed. For example, one may assume that the mixture is composed of two 
Gaussians even though their means and standard deviations are unknowm. To es- 
timate the means and the standard deviations, additional assumptions will have 
to be made before data can be utilized for obtaining the Gaussians. A “K-means” 
algorithm described in Reference [37] can be used for this purpose. 

Of the two broad categories of clustering methods, unsupervised clustering 
is more useful in practice. This is due to the following factors: (1) for certain 
problems it is not easy to label the training samples due to their size, (2) the 
clusters can undergo small changes, and (3) very little is often known about the 
structure of the data. One of the ways of discovering structures in the data is 
by constructing a weighted graph. Distance relationships in the graph can then 
be used to partition the graph into sub-graphs to further improve the distance 
relationships. The graph-theoretical method described in Reference [104] uses a 
minimum spanning tree to partition the set of points into perceptually organized 
clusters. The perceptual organization is defined by the principles of proximity, 
similarity and continuity. 

As an example of unsupervised clustering, Reference [95] describes a hierar- 
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chical clustering method for grouping discrete ranges for the scene shown in Figure 

2.10. The technique described there first represents the range histogram as a sum 
of Gaussian. Next, the features are grouped based on separation in the horizontal 
plane. Finally, an algorithm based on the minimum spanning tree (MST) [104] 
is used for grouping the range points based on the separation in the image plane. 
The results from unsupervised clustering for a sample scene is shown in Figure 

2.11. It can be observed from the figure that the method generated six clusters. 



Figure 2.11: Groups in the image using unsupervised clustering. 


Another way of addressing the clustering problem is to cast it as a discrete 
optimization problem which minimizes a certain distance function. Distance func- 
tions such as within-cluster and between-cluster distance measures based on the 
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scatter matrices are discussed in [27]. Since the set of features or points is finite, 
there can only be a finite number of partitions. Thus, in theory, the clustering 
problem can always be solved by exhaustive search. However, in practice, such an 
approach is not feasible because there are approximately c n /c! ways of partitioning 
a set of n elements into c groups. Due to this reason, the approach most frequently 
used is that of iterative optimization. 

In Reference [14] application of the Monte-carlo methods for clustering range 
points into objects is described. It should be noted that these methods guaran- 
tee local but not global optimization [27]. Despite these limitations, the fact that 
computational requirements are reasonable make these approaches desirable. A 
technique based on Simulated Annealing for refining the initial grouping is de- 
scribed in [51]. The initial grouping in this case is obtained by assigning the range 
points to image regions obtained by labeling a segmented image. 

2*3.5 Geometric Modeling 

The clustering of discrete range points enables one to assume the range to 
be continuous within a group. It is possible to subsequently create a dense range 
map via interpolation within the groups. Modeling of dense range images has been 
studied by several authors [10, 58, 61]. The dense range images can be modeled 
into objects by fitting surfaces using polynomials, splines [74], Delaunay triangles 
[23] and other mathematical surfaces. 

Several different approaches for representing surfaces defined by a set of ran- 
domly located points using triangular grids are described in [23]. These representa- 
tions approximate the surface as a network of connected triangles with vertices at 
the data points. Many of the surface fitting algorithms use the properties of Delau- 
nay triangles to discretize the domain with triangular elements. These algorithms 
may be broadly classified as incremental algorithms and divide-and-conquer al- 
gorithms. Incremental algorithms start from a boundary or interior point and 
create triangles by adding the remaining points. Divide-and-conquer algorithms 
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recursively split the set of data points into equal subsets until elementary sets are 
obtained, and then merge them pairwise. For example, application of the incre- 
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Figure 2.12: Surface representation of the groups in Figure 2.11 by triangular 
elements. 

mental algorithm in Reference [64] to the clusters in Figure 2.11 yields a surface 
representation shown in Figure 2.12. 

Once object modeling is accomplished by surface representation, additional 
geometric details can be extracted using surface interpolation. Several different 
element shapes and shape functions are discussed in the Finite Element Method 
literature [80, 105]. Some of these can be used for efficient interpolation. For the 
example shown in Figure 2.12, the interpolated range data is encoded as gray- 
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levels and presented in Figure 2.13. The figure also shows the modeled ground 



Figure 2.13: Finite element object models for surfaces in Figure 2.12. 


plane. 

For ground plane modeling, a Least Squares method can be used with the 
points from every group that are below a certain altitude. An assumption implicit 
in such modeling is that all the objects observed in the scene lie on a ground plane. 
A perspective projection of a rectangular grid on the ground plane can then be 
created to aid visualization. This process requires knowledge of camera altitude, 
and pitch and roll angles with respect to a local horizontal. Using the relative 
geometry of the camera with respect to the ground, the locations of horizon and 
the vanishing point can be obtained. The relationships needed for obtaining the 
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ground plane representation are described in [27] and [44]. The representation of 
ground plane using a grid projection is shown in Figure 2.14. The locations of the 
horizon and the vanishing point are also shown in the figure. The grid size of 12.5 



Figure 2.14: Ground plane representation, 
feet by 20 feet was used in this case. 

For the example presented above, inference is direct once the scene is assumed 
to be a model of a plane with objects lying on it. Geometric modeling directly 
yields the orientation of the plane with respect to the local horizontal, and the size, 
distance and shape of the objects with respect to the camera. Note that general 
scene understanding is much more complex and requires sophisticated inference 
techniques. 


47 


2.3.6 Inference Techniques 

In a vision system, once the primitives or features are derived and classified, 
rules of inference can be used for object recognition. Many vision systems implicitly 
assume object models in order to aid the object recognition process and to develop 
an understanding of the scene. Comparative studies of several model-based object 
recognition algorithms are discussed in References [12] and [20]. 

Current model— based object recognition systems have several limitations. One 
of them is the difficulty in representing and describing objects. Only simple objects 
can be recognized by matching two-dimensional features with two-dimensional 
object models. The non-availability of higher dimensional features restrict the 
recognition capabilities to few object classes viewed in a particular way. A more 
general system will require the ability to extract three-dimensional features that 
are view point independent and match them with three-dimensional object models. 
Another difficulty is the non-availability of descriptors of surface properties of 
objects. Three issues that a model-based object recognition system has to deal 
with are: (1) design of features that describe physical properties and their spatial 
relationships, (2) a meaningful representation of the feature vector for an object 
class and (3) matching between the feature vector and object models for object 
recognition in a general scene [20]. 

The discussion of model— based object recognition with the background of 
algorithms and processes discussed in the previous sections indicates that it is 
difficult to design a general purpose vision-based object recognition system, and 
that a sequence of several low-level and high-level vision techniques are needed. 
Finally, without a model, the task of object recognition is virtually hopeless. In 
view of these observations, model— based techniques appear to offer the most direct 
scene interpretation without using elaborate inference techniques. 
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2.3.7 Ranging 

Recovery of three-dimensional geometry from two-dimensional images is based 
on the fact that the differences between the locations of objects in two or more 
images obtained from different vantage points is a measure of their range. The pro- 
cess of finding the same object in multiple images is known as the correspondence 
problem in the machine vision literature. Since real images have limited field-of- 
view and resolution, the concept of correspondence is effectively a correspondence 
hypothesis. The relative object displacement obtained by satisfying this hypoth- 
esis in the image plane is called disparity. Due to perspective projection from 
the three-dimensional scene to the two-dimensional image, the farther the object 
is from the imaging device, the less disparity it exhibits. Closer objects exhibit 
larger disparity. Many vision— based methods discussed in the literature compute 
the disparity, thus recover the range to objects in the scene [7, 11, 46, 66, 81, 90]. 

In the simplest case of stereo vision where a pair of images are acquired by 
two cameras separated by a baseline, range can be computed by triangulation. 
For example, consider the geometry in Figure 2.15. In this figure, a point object 
appears along the line connecting the camera centers at r/i in one image and at 
u 2 in the other image. Let the distance between the camera centers be b and the 
camera focal length be /. The Azimuth angles with respect to the optic axes of 
the two cameras then are ?/’i = tan~ l (u\/ f) and i[> 2 = tan *( 1 ^ 2 / /)• Since the two 
angles and the base of the triangle are known, the lengths of the range with respect 
to the cameras can be computed. An equivalent calculation can be done using a 
single moving camera. In this case, motion establishes the baseline required for 
triangulation. This is also known as Cyclopean vision, inspired by the mythical 
single eyed monster in Homer’s Odyssey. 

Driven by the needs of helicopter nap-of-the-earth guidance problem, ma- 
chine vision techniques for ranging has recently attracted significant research at- 
tention. There are two distinct classes of algorithms that determine range b^ 
satisfying the correspondence hypothesis. They are known as field-based and 
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Figure 2.15: Range determination by triangulation. 
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feature-based methods. Field-based methods, such as [7, 46, 66, 81], assume a 
continuous variation of image intensity as a function of position or position and 
time. Feature-based methods, such as [18, 60, 76, 90], identify features in images, 
such as points, lines and contours in order to establish correspondence. Both, 
field-based and feature— based methods can only compute range to few locations 
in the image. Although field-based schemes have the potential of providing denser 
range maps, experience has shown [67, 68] that range can be reliably computed 
only at about 10% of the points in the image. This is due to the fact that compu- 
tations break down in regions of near uniform brightness. Feature-based methods 
by their very nature can only compute range at discrete locations. An example of 
the range computations using a feature-based method is given in Figure 2.10. 

Reference [66] describes a field-based ranging procedure using motion se- 
quences generated by a single camera fixed to a moving vehicle. The method 
is based on the Optical Flow Constraint Equation of Reference [46] that relates 
the temporal partial derivatives with the spatial derivatives of the image function. 
Due to the use of partial derivatives, a smoothness constraint has to be enforced for 
the computation of range [46]. As discussed in Reference [66], incremental perspec- 
tive projection equations can be directly combined with Optical Flow Constraint 
to yield a single navigation equation. This equation can then be used for obtaining 
the range. A crucial part of this method is evaluation of partial derivatives of the 
image function. In [66], the partial derivatives are estimated using a method based 
on the Calculus of Variations. 

In Reference [67], a ranging scheme using image pairs is described. A multi- 
dimensional Taylor Series approximation of the correspondence hypothesis is used. 
The advantage of this method is that it does not require temporal partial deriva- 
tive of the image function. Hence, this formulation does not require the concept 
of optical flow. In this method also, the incremental perspective projection equa- 
tions are used with the Taylor series approximations to formulate a navigation 
equation. Since the navigation equation depends on the spatial partial derivatives, 
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the derivatives are computed using a finite difference scheme. Camera transla- 
tion distance between the image pairs is used to perform range computations. No 
rotational motion is assumed. 

To overcome the difficulties associated with derivative computation of noisy 
image functions using finite difference schemes, causal estimators that attenuate 
the noise in the process of derivative estimation are proposed in Reference [68]. 
This method offers the possibility of derivative computation during the image 
data collection process. Multi-dimensional Taylor series approximation of the 
correspondence hypothesis is used in this study also. Eliminating the disparities 
in favour of camera motion parameters and scene depth using the perspective 
projection equations, an Optical Ranging Polynomial is obtained. This polynomial 
is then solved to obtain the range. The algorithm has been demonstrated on a 

stereo image pair of a laboratory scene. 

An extension of the previous algorithm that includes both translational and 
rotational displacements is described in Reference [70]. The central theme, in- 
cluding the technique for estimation of the partial derivatives, is the same as the 
previous algorithm. The algorithm has been demonstrated on an outdoor image 
sequence acquired by a camera mounted on the nose of a helicopter. The im- 
ages are temporally separated and were acquired as the helicopter underwent both 
translational and rotational motion. 

To overcome the difficulties of noise amplifying derivative estimation process, 
a derivative free ranging method is proposed in Reference [71]. In this algorithm, 
the correspondence hypothesis is approximated using Pade’ approximation and 
used as a differential constraint in an optimization problem with a quadratic cost. 
The state variable is the sum of image functions of the two stereo images and the 
control variable is the range to objects seen in the images. The resulting necessary 
conditions for optimality are linear permitting the solution using the backward 
sweep method [13]. The method was demonstrated on a pair of stereo images of a 


laboratory scene. 
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Feature-based ranging methods are considered next. Given corresponding 
features in two images, the study in [90] describes how the estimated object location 
is influenced by the estimation algorithm and the relative geometry between the 
camera and the object. Three different Kalman Filter formulations are proposed 
for ranging in Reference [90]. These are: (1) inertial coordinate formulation (2) 
sensor coordinate formulation and (3) polynomial model for image point motion. 
The polynomial filter was found to be unsuitable for general camera motion. The 
methods described in Reference [90] assume that the motion parameters such as 
camera position, attitude, and translational and rotational velocity are available 
from an onboard Inertial Navigation System. 

Research reported in Reference [86] develops a normalized correlation function 
based feature correspondence procedure. This technique forms the first step in 
the Kalman filtering algorithm. Feature detection is accomplished by using an 
edge operation and correspondence is achieved by using the gray-levels of the 
detected features. In addition, a recursive algorithm for range estimation based 
on translational motion is also described. Since translational motion is assumed, 
the search for correspondence is restricted to envelopes along radial lines eminating 
from the focus-of-expansion in the images. Results for a laboratory image sequence 
are obtained by using the recursive algorithm. 

Details of the correspondence procedure when images are acquired from a 
camera undergoing general motion are described in Reference [88]. Thus, this work 
extends the procedure given in Reference [86] for more general motion involving 
translation and rotation. An elliptical search window based on the propagated 
range estimate is used to minimize the search effort. Results for a laboratory 
image sequence are described. 

Results for an outdoor image sequence obtained with the Kalman Filter for- 
mulated in the sensor frame in Reference [90] are reported in Reference [89]. It is 
shown that good range accuracy is obtained for the objects in the field-of-view of 
the camera. This result is significant because it is much more difficult to establish 
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correspondence in the images taken from a rotating and translating platform, such 
as a helicopter in flight. Details of the image acquisition procedure are described 
in Reference [82]. The results of application of the motion algorithms on the lab- 
oratory and flight image sequences are described in [91]. That report summarizes 
the procedure and the results of References [88] and [89]. 

Since the objects in the field-of-view are at various ranges, It may be ad- 
vantageous to use a different measurement rate in different portions of the image. 
This idea is explored in Reference [92]. The technique for range estimation in- 
volves accepting the measurement for the Kalman Filter only when the tracked 
feature moves more than a set threshold in the next image. Numerical results 
presented for the outdoor image sequence show that the multirate filter provides 
the same estimation accuracy as the standard Kalman Filter, with a significantly 
lower computational effort. Since different features are updated at different times, 
the book-keeping task is more involved when compared to the single rate filter 
implementation. 

When range information is obtained using a single camera, it is sensitive to 
the direction of motion. Hence, the estimates are poor close to the focus-of- 
expansion. An analysis of motion and stereo methods is provided in Reference [87] 
to demonstrate that motion methods provide more accurate range information 
away from the focus-of-expansion and stereo methods provide superior accuracy 
close to the focus-of-expansion. In order to overcome the limitations of the stereo 
method, a recursive stereo method is described in Reference [87]. This method 
is then contrasted with standard stereo method and the earlier recursive motion 
algorithm [86]. It is suggested that an integrated stereo and motion method based 
on the recursive motion method and the recursive stereo method has the potential 
for providing more accurate range estimates when compared to either of the two 
methods. 

A hybrid motion/stereo algorithm is described in Reference [83]. This algo- 
rithm is an extension of the recursive motion algorithm given in Reference [90]. 
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The range predictions generated by the Kalman filter are used for constraining the 
search space for feature correspondence in the stereo and motion pairs. One of 
the advantages of the method is that the Kalman Filter can be initialized with the 
range estimates obtained by processing the stereo pair. Results of application to 
an outdoor image sequence shows that the hybrid estimates are an improvement 
over the monocular estimates. Both methods generate estimates which appear to 
converge to the true range over time. 

These algorithms have been applied to numerous images of outdoor scenes ob- 
tained from helicopter-borne cameras. The results obtained using these algorithms 
have been verified against range data obtained by a laser range finder. These al- 
gorithms can be considered to represent a mature class of vision based ranging 
algorithms. 

The following conclusions can be drawn based on the review of the field-based 
and feature-based methods: (1) correspondence of regions in one image to regions 
in another image is the most significant problem, (2) even in an unstructured 
scene, ranging algorithms can be made to work reasonably if the camera position 
and angular displacements are available from an independent source like an Inertial 
Navigation System, (3) inclusion of system dynamics in the design of a recursive 
state estimator leads to higher estimation accuracy, and (4) a hybrid motion- 
stereo method provides higher accuracy when compared to pure motion and stereo 
methods. 

Egomotion or self motion problem is the dual of the ranging problem. In 
this case, the camera position and orientation are the unknowns to be determined. 
In order to solve the problem, it is often assumed that the objects in the field— 
of-view are stationary. Along with this assumption, if correspondence can be 
established between features in successive images, the change in camera position 
and orientation can be computed. To determine the absolute camera position and 
orientation, the location of objects in the field-of-view have to be known with 
respect to an inertial coordinate system. Thus, the solution of egomotion problem 
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requires an underlying scene model via the definition of an inertial coordinate 
system together with the locations of the objects with respect to this coordinate 
system. This suggests two possible approaches: (1) use of the vision based ranging 
algorithms to recover the absolute camera position and orientation, or alternatively 
(2) use of a model-based approach to directly recover the camera position and 
orientation. 

Camera calibration problem [6, 27] is a special model-based approach in which 
the correspondence between the objects in the scene and the image are known. The 
objective of camera calibration is the determination of camera optical character- 
istics and the camera position and orientation. Intrinsic camera parameters are: 
focal length, lens distortion, scale factor, and center of the image plane. In order 
to calibrate the camera, a planar grid target is placed at a certain orientation and 
distance away from the camera, and an image of the grid is obtained. The position 
of every grid point on the plane is known with reference to an inertial frame. In 
addition, the correspondence between every grid point on the plane and in the 
image are known. Since every grid location of the plane is related to its image via 
the intrinsic and extrinsic parameters, these parameters can be computed using 
an iterative algorithm. The traditional approach is to use a nonlinear parameter 
optimization technique. An alternative two-stage technique for camera calibra- 
tion is described in Reference [98]. This technique solves the problem by using the 
least-squares method. Only few parameters are computed using nonlinear search. 
Efficiency of the process can be greatly enhanced by generating the initial guess 
using the least-squares method. 

If correspondence is unknown, the calibration technique is not applicable even 
with known intrinsic parameters. Model-based methods are useful in this case. A 
model-based method that uses local feature correspondence and a Kalman Filter 
is described in References [25] and [26]. In References [25] and [26], the initial 
camera position and orientation estimates are used together with the perspective 
projection equations for projecting known model features such as corners and curve 


segments into the image plane. This results in the creation of the model image. 
These features are then identified in the actual image acquired by the camera. The 
position difference of the features in the actual image and the model image are used 
to drive the Kalman Filter to improve the camera position and orientation estimate. 
Locating model features in the actual image is the main limitation of the algorithm. 
The difficulty is caused by the fact that in general, the features in the actual image 
are significantly different than the features in the model image. Therefore, the 
search involves matching a considerably simplified model template with regions in 
the image. In contrast, the ranging methods reviewed earlier use a templete based 
on a previous real image. In order to work satisfactorily, model based matching 
requires the features to be invarient to scale and rotation. It is difficult to find such 
features in real scenes. If the scene is such that one feature cannot be distinguished 
from another, matching the features may be difficult because a model feature could 
potentially match with many image features. 

2.3.7. 1 Uniqueness of Solutions 

Before embarking on the development of pilot aids based on machine vision 
techniques examined in this chapter, it is essential to address the question of 
uniqueness of the solutions. Machine vision literature [48, 73] poses the uniqueness 
problem as follows: given displacements and velocities of image points, under what 
conditions is it possible to recover the shape of the scene and the relative motion 

between the camera and the objects in the scene? 

For differential motion, the research given in Reference [48] shows that ambi- 
guities arise only in the case of certain hyperboloids of one sheet and their degen- 
eracies, such as circular cylinders, elliptic cones, hyperbolic paraboloids, and two 
intersecting planes that are viewed from a point on their surface. The governing 
equation of hyperboloids of one sheet is [85]: 
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where (x, ?/, z) is the coordinate of a point on the surface of the hyperboloid, and a, 
b and c are constants. For large motion, it is shown in Reference [28] that if only five 
points are available in the image, then up to ten solutions are possible. Research 
in Reference [97] shows that the solution is usually unique if the displacements of 
seven points in two successive images are known. The solution is non-unique only 
if these points lie on a cone passing through the origin or on two planes, with one 
plane passing through the origin. 

More recently, it has been shown in [73] that, only certain hyperboloids of one 
sheet and their degeneracies when viewed from a point on their surface can give 
rise to ambiguity. Moreover, Reference [73] shows that in the case of hyperboloids 
of one sheet and hyperbolic paraboloids, there can be at most three solutions. 
That work also demonstrates that in the case of intersecting planes and circular 
cylinders, there can exist at most two solutions. It is also pointed out that cones 
cannot give rise to ambiguity unless the motion is differential. 

The next chapter will apply the machine vision techniques described in this 
chapter to develop pilot aids for night landing. Machine vision algorithms for air- 
craft position, velocity and attitudes with respect to the runway will be derived 
and evaluated in simulations. Data sources for these algorithms will also be iden- 
tified. Past research examined in this chapter indicates that it may be possible 
to synthesize machine vision systems that produce unique solutions for pilot aid- 
ing during night landing. This conclusion is a consequence of the two facts: (1) 
the airport lighting layout is viewed from above and (2) the underlying lighting 
geometry is planar. 

2.4 Summary 

Since the algorithms for machine vision are subject areas of Computer Vision, 
a review of the literature relevant for the design for such algorithms was discussed 
in this chapter. Two broad classes of algorithms, low-level and higher-level vision 
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algorithms were examined. Several algorithms for low-level vision process of im- 
age conditioning and edge detection were outlined. Higher— level vision algorithms 
for boundary detection, gray-level segmentation, texture segmentation, clustering, 
geometric modeling, inference and ranging were also discussed. Results of applying 
several of the low-level and high-level algorithms to an actual image of a runway 
scene were presented. These results illustrate the types of information that can 
be derived using image processing algorithms. Since position determination is the 
central topic of this research, field-based and feature-based algorithms reported in 
the literature that are closely related to this work were described. Many of these 
methods required the position and orientation of the camera to be known. As a 
result, they are not directly useful for the runway relative aircraft position and 
orientation determination problem. The literature for camera calibration problem 
was also found to be inadequate because those techniques assume correspondence 
between the grid points on the calibration plane and the image plane. Since corre- 
spondence between the model of the runway scene and the image acquired by the 
camera is unknown a priori, template— based local feature correspondence methods 
were also found to be unsuitable. Finally, the question of uniqueness of solutions 
was addressed. Based on the available literature, it was established that a unique 
solution of the runway relative position and orientation could be found for the 
viewing geometry used in this research. 



Chapter 3 


Machine Vision Based Landing 
Aids 


This chapter develops the basic building blocks for constructing the machine 
vision algorithms for aircraft runway relative position and orientation estimation. 
With this goal, the nature of the landing task and the accuracy requirements are ex- 
amined first. Clearly, the specification of a runway fixed inertial coordinate system, 
body coordinate system and camera coordinate system are essential components 
of methods for aircraft position and orientation estimation. Aircraft equations of 
motion which relate the aerodynamic and propulsive forces and moments to the 
translational and rotational motion of the aircraft are then developed with respect 
to the body and inertial coordinate systems. An onboard pinhole camera model 
that relates the inertial location of an airport light to its location in the image 
plane is described subsequently. Finally, landing and image simulation procedures 
are described to tie these building blocks together. 


3.1 Aircraft Landing Operation 

Aircraft arrival flight to the destination airport can be broken up into two 
broad segments: en-route descent and final approach to touchdown. A host of 
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procedures have to be followed by the pilot during both these phases. During 
descent from cruise, the rate of descent and airspeed have to be controlled to 
comply with the restrictions imposed by the air traffic control system. The airspeed 
has to be reduced to 250 Knots Indicated- Airspeed or less, when below 10,000 feet 
MSL (Mean Sea Level) [24]. 

A published arrival procedure called Standard Terminal Arrival (STAR) is 
used to transition from the en-route structure to an outer fix or an instrument 
approach fix or an arrival waypoint in the terminal area. The aircraft is then flown 
to the final approach fix to intercept the glide slope. The optimum length of the 
final approach is five miles; the maximum length is ten miles [24], Once on the glide 
slope, aircraft speed, rate of descent and certain altitude distance relationships are 
maintained until the aircraft is beyond the runway threshold and at a prescribed 
altitude. At this stage, the aircraft executes the flare maneuver to achieve a gentle 
touchdown. 

During landing, the pilot controls the aircraft lateral displacement from the 
runway centerline, distance from the touchdown point, altitude, yaw-pitch-roll 
orientations, rate of descent and rate of closure with the touchdown point on the 
runway surface. The desired glide path which describes the altitude, time and 
distance relationships during a typical landing are shown in Figure 3.1. Figure 3.1 
shows a commonly employed three degree glide slope approach. Glide slope, r*, is 
defined as: 

v = tan ^(h/xgo) (3-1) 

where h is the altitude and x go is the distance-to-go to the touchdown point. 
For precision approach, the glide slope is between 2.5 and three degrees at most 
airports [24]. The glide slope in conjunction with the location of the touchdown 
point specifies the desired aircraft position with respect to the runway threshold 
as a function of altitude. Optimal threshold crossing height is 50 feet but it may 
be as high as 60 feet or as low as 32 feet [24]. The touchdown point is specified in 
terms of the distance from the runway threshold. From Figure 3.1, it may be seen 
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Figure 3.1: Glide slope, altitude, time and distance relationships. 

that three degree glide slope requires the aircraft to be at distances of 6633 feet, 
2816 feet, and 908 feet corresponding to altitudes of 400 feet, 200 feet and 100 feet 
respectively. These distances translate to 30 seconds, 13 seconds and four seconds 
to the runway threshold. Time-to-go calculations are based on a typical approach 
speed of 220 feet /second (130 knots). 

The discussion of landing procedures is incomplete without mention of the 
abort procedures. Once the aircraft has passed the final approach fix, it is flown 
to the minimum descent altitude with enough time and distance remaining to 
identify the runway environment before continuing on the visual approach to the 
touchdown point. Descent below the minimum descent altitude is not authorized 
until visual reference with the runway environment is established and the aircraft 
is in a position to execute safe landing [24], If it is unable to execute a safe 
landing, the aircraft is flown at or above the minimum descent altitude to the 
missed approach point. Subsequently, the aircraft is routed back to the outer fix 
for another landing attempt. Depending on the ground and airborne equipment, 
the decision to land can be delayed as the aircraft is flown along the glide slope. 
There are prescribed landing categories with associated decision heights up to 
which aircraft can be flown with instruments. Beyond the decision height for 
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the category of landing, it should be possible to fly the aircraft with just the 
visual reference. Like the final approach fix data, missed approach points are also 
published in navigation charts. 

3.1.1 Landing Accuracy Requirements 

One of the terms used for landing conditions is runway visual range (RVR). 
RVR is the distance from which the pilot can see the high-intensity runway edge 
lights. It is determined by transmissometer measurements near the threshold [49]. 
The transmissometer consists of a light source with a narrow beam projector and a 
receiver with a narrow beam acceptance angle. In order to make the measurements, 
these two components are raised to 15 feet above ground and separated by 500 
feet. The amount of light received is a measure of atmospheric transmissivity. The 
measurements are compensated for the intensity setting of the edge lights and the 
time of day or night. For category II and III operations, two measurements are 
made. One near the threshold and the other near the midpoint of the runway. 
While useful, these measurements do not accurately predict the visibility along 
the approach path, since the measurements are made close to the ground. 

Visibility on the runway is classified into I, II and III categories. Category III 
is further subdivided into a, b, c. The three categories are defined in terms of the 
RVR and decision height. Decision height (DH) is defined as the minimum height 
above the runway where a decision must be made by the pilot to continue descent 
to landing or to abort. The decision is based on the pilot being able to obtain 
visual guidance cues provided by airport lighting without depending on cockpit 
instruments. 

The various categories and the associated RVR and decision heights are listed 
in Table 3.1 [31, 32]. Capability for automatic landing all the way to touchdown 
is required for all category III landings. For category Ilia, the rollout after land- 
ing and taxiing is manual. For category Illb, an automatic rollout capability is 
additionally required. For category IIIc, an automatic taxiing capability is also 
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Table 3.1: Visibility Categories 


Category 

Decision Height (ft) 

Visibility (ft) 

CAT I 

200 

1800 

CAT II 

100 

1200 

CAT Ilia 

50 

700 

CAT Illb 

0 < DH < 50 

150 

CAT IIIc 

0 

0 


required, in addition to the landing and rollout capabilities. 

Although Table 3.1 lists the RVR and decision heights for the various cat- 
egories, it does not list the navigation accuracy requirements. The performance 
specifications for Federal Aviation Administration (FAA) defined precision ap- 
proach and landing categories are given in Table 3.2 [100]. These accuracy re- 

Table 3.2: Aviation Navigation Accuracy Requirements 


Category 

Lateral (ft) 

Vertical (ft) 

CAT I 

±56.1 

±13.45 

CAT II 

±17.06 

±5.58 

CAT III 

±13.45 

±1.97 


quirements will be used to evaluate the performance of the machine vision based 
algorithms developed in this report. 


3.2 Coordinate Systems 

Various coordinate systems used in this report are illustrated in Figure 3.2. 
In this figure, i is the origin of the inertial frame attached to the runway threshold. 
Since the location of all lights are given with respect to the threshold bar, it 
is a natural choice for the location of the origin of the inertial coordinate system. 
Furthermore, since the centerline lights form a principal axis of symmet ry, following 
the flight dynamics convention, the x-axis of the inertial frame is aligned with the 
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runway centerline in the approach direction and the z-axis points down. The y- 
axis completes the right-handed triad. The origin of the aircraft body axes is 
located at the point b. Its position with respect to the inertial frame is given 
by the vector X* b with components x b , y b and z b . The camera frame is located 
at c. The camera position with respect to the body frame is given by the vector 
Xj? with components l x , l y and l z . Since the camera is rigidly attached to the 
aircraft structure, the vector X]? is assumed to be constant in the present research. 
Let p be a light on the runway and let its position with respect to the inertial 
frame be given by the vector Xp with components x p , y v and z p . Also, let the 
position of point p with respect to the camera frame be given by the vector Xp 
with components x cp , y Cp and z cp . 

The position of the point p with respect to the aircraft in the inertial frame 
is given by the vector 

X^ - X[> = [(x p - Xi,), (y p - y b ), (z p - z b )] T (3.2) 

The transformation matrix from the inertial frame to the body frame Tb/i can be 
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obtained in terms of the yaw attitude pitch attitude 0, and roll attitude <f) as 
[96]: 


Tb/i = 


cos ip cos 6 sin ip cos 0 — sin 6 

— sin ip cos (p + cos ip sin 0 sin <j> cos ip cos <j> + sin ip sin 0 sin <p cos 8 sin <f> 

sin ip sin <p + cos ip sin 6 cos <p - cos ip sin <f> + sin ip sin 0 cos <p cos 6 cos <p 

(3-3) 


The position of point p with respect to the body frame can be obtained as: 



T b/i (Xt - Xt) 


(3-4) 


Similarly, the position of point p with respect to the camera frame is given by: 

X' = T c/b (Xp - Xj) (3.5) 

Here, T c /b is the constant transformation matrix from the body frame to the 
camera frame. Combining Equations (3.4) and (3.5): 

X' = T c / b T b /i(Xp - Xj.) - T c / b X b (3.6) 

Since the camera is assumed to be fixed with respect to the body, the product: 
— T c /bX|? is a known constant vector k with components k x , k y and k z . Further- 
more, if through rg are defined as the elements of the transformation matrix 
from the inertial frame to the camera frame, T c/ j = T c /bT b /i, the components of 
the position vector Xp can be obtained as: 

x cp = ri(x p - Xb) + r 2 (y p - Vb) + r 3 (z p - z b ) + k x (3.7) 

y Cp = r 4 (z p - Xi,) + r 5 (y p - 7/t ) + r 6 (z p - z b ) + k y (3.8) 

z cp = r 7 (x p - x b ) + r 8 (t/ p - y b ) + r 9 (z p - z b ) + k z (3.9) 

Equations (3.7) through (3.9) show the relationships between the location of 

the airport lights in the camera frame and the aircraft position and orientation. 
The aircraft position and orientation evolve due the to forces and moments acting 


on it. 
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3.3 Aircraft Dynamic Model 

Aircraft is subjected to aerodynamic, gravitational and propulsive forces and 
moments. These forces and moments result in translational and rotational motion 
of the aircraft. Equations of motion relate the airplane motion to the forces and 
moments. Three coordinate systems are used to express the forces and moments 
in a convenient way. These are described first. Subsequently, equations of motion 
are presented with forces and moments expressed in these coordinate systems. 


3.3.1 Coordinate Systems 


Aerodynamic forces and moments on the aircraft depend on the orientation 
of the airframe with respect to the airflow. Since rotation around the free-stream 
velocity vector in a uniform airflow does not cause changes in the aerodynamic 
forces and moments, they depend only on two orientation angles with respect to 
the relative wind. These are the angle of attack, a and the angle of sideslip, 0 
illustrated in Figure 3.3. 

Figure 3.3 shows the body axes system with the x-axis aligned with the fuse- 
lage reference line, the z-axis in the aircraft plane of symmetry and orthogonal to 
the x-axis, and the y-axis normal to the plane of symmetry. The angle of attack 
and sideslip are defined by performing plane rotation about the body y-axis by 
a, followed by another plane rotation about the new z-axis by 0 such that the 
x-axis is aligned with the relative wind. The variables a and 0 are the angle of 
attack and angle of side slip respectively. The axis system resulting from the first 
rotation about the y— axis is often called the stability axis system. 

With the angles of attack and sideslip defined by the axes systems, the trans- 
formation from the body to stability axes T s/b is given by: 


Ts/b — 


cos a 0 sin a 
0 1 0 
— sin a 0 cos a 


(3.10) 
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Figure 3.3: Aircraft axes and angles. 

The transformation from the stability to the wind axes T w / 8 is given by: 

cos 0 sin 0 0 

-w/s = — sin 0 cos 0 0 (3-11) 

0 0 1 

Concatenating the two transformations, the transformation from the body axes to 
the wind axes T w /b can be obtained as: 

cos a cos 0 sin 0 sin a cos 0 
— cos a sin 0 cos 0 —sin a sin 0 (3.12) 

— sin a 0 cos a 

Forces and moments can be expressed in body or wind axis systems by using 
these transformation matrices. 


- w/b 


= T w / 8 T s/ b = 


3.3.2 Forces and Moments 

The forces and moments on the aircraft arise due to aerodynamics, gravita- 
tional acceleration and engine. The aerodynamic forces are specified in the wind 
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axes. The components along the negative x-axis, positive y-axis and the negative 
z-axis are called drag, Z), sideforce, Y and lift, L. In order to avoid dealing with 
state dependent moments of inertia, the aerodynamic moments are defined in the 
body axes. These are: the rolling moment L r about the x-axes, the pitching mo- 
ment M about the y-axis, and the yawing moment N about the z-axis. The forces 
and moments are specified in terms of dimensionless coefficients as follows. 


D = qSC D 

(3.13) 

L = qSC L 

(3.14) 

Y = qSCy 

(3.15) 

Lr = qSb s C Lr 

(3.16) 

M = qSb c C M 

(3.17) 

N = qSb s C N 

(3.18) 

where, q is the free-stream dynamic pressure, S is the 

wing reference area, b s is 


the wing span, b c is the wing mean aerodynamic chord, Co is the drag coefficient, 
Cl is the lift coefficient, Cy is the sideforce coefficient, Cl t is the rolling moment 
coefficient, Cm is the pitching moment coefficient, and Cn is the yawing moment 
coefficient. The aerodynamic coefficients primarily depend on the aerodynamic 
angles, a and 0, the Mach number and control surface deflections. They also 
depend upon the body rates. A detailed discussion of these coefficients is available 
in Reference [96]. The forces and moments due to the engine arise from the thrust, 
its location with respect to the aircraft center of gravity and misalignment angles. 
Thrust-related forces and moments are denoted by subscript T in the following. 

Forces and moments due to aerodynamics and engine thrust in the body axes 
are: _ _ 


F x 


-D 


i 

£ 

i 

Fy 

_ rp- 1 

x w/b 

Y 

+ 

Fy T 

F z 


L 


Fz t 


T b = 


(3.19) 
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and 



u 


Lr 


Lrj' 

TB = 

Mb 

— 

M 

+ 

Mj 


. Nb . 


N 


Nj 


(3.20) 


Here, F B is the force vector in the body axes, F Xt , F Yt and F Zt are the components 
of the engine thrust vector in the body axes. T w /b i s the transformation matrix 
from the body to the wind axes, described earlier in Equation (3.12), and t b is the 


moment vector in the body axes. The variables T rr , Mx and Nx are the rolling 
moment, pitching moment and yawing moment due to the engine thrust. 


3.3.3 Equations of Motion 

Six-degree-of-freedom aircraft model is given by twelve first-order nonlinear 
differential equations involving the position and attitude dynamics. Ihe twelve 
state variables consist of: (1) the inertial position of the aircraft represented by its 
topocentric coordinates x,y,z] (2) aircraft velocity components measured in the 
body axes U, V, W ; (3) the body Euler angles denoted by in roll, pitch and 

yaw, respectively; and (4) the angular rate of the body p, in body axes. The 
aircraft equations of motion are given in the following. 

The force equations are [96]: 


„ Fx 

U = rV — qW — q sin 8 4 

(3.21) 

m 

V = — rU + pW + q sin 6 cos 0 H — - 

(3.22) 

m 

W = qU — pV + g cos <f> cos 0 4 - 

m 

(3.23) 


Here, g is acceleration due to gravity and m is the mass of the aircraft. 
The rotational kinematic equations are [69]: 


t/i = q sin 4> sec 8 + r cos <f) sec 6 

8 = q cos <f> — r sin <f> 

cf) = p -f q sin <j> tan 9 + r cos <j>t&n 9 


(3.24) 

(3.25) 

(3.26) 
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The angular acceleration equations are [69]: 

p = \I\Lk 4 7 2 M fc 4- IjNb — p 2 {hzh ~ hyh) 

-|- pq(I X zh ~ hzh ~ hh) - pr{h y h + hh ~ hzh) 

4 q 2 (I y zh — I xy h) ~ qr{hh — hyh + hzh) 

- r 2 (l yz I\ - hzh))/det(I) (3.27) 

q = [ /'i I-b 4 hMb 4 hNb — p 2 (I rz 1 4 h: y h ) 

4 pq(hzh ~~ lyzl 4 — ^ 9 / 5 ) — P r(I xy h 4 I&I 4 hzh) 

4 q 2 (I y zh ~ hyh) ~ qr{hh ~ hyh + hzh) 

- r 2 (I y J 2 -I X zh)]/det(I) (3.28) 

r = [7 3 7, 6 4 hMb 4 hNb ~ p 2 {hzh ~ hyh) 

4 pq{hzh ~ lyzh ~ hh) ~ pr(h y h + hh ~ hzh) 

4 q 2 {Iyzh ~ hyh) ~ qr(hh hyh T hzh) 

- r 2 (Iyzh - Ixzh)]/det(I) (3.29) 

where the determinant of the inertia matrix is given by: 

det(I) = I x Iyh - Zhyhzlyz ~ hh ^ ~ Iyhz 2 ~ hhy (3-30) 

and 

7! = Iyh — lyz 2 ( 3 - 31 ) 

h = hyh + Iyzhz (3-32) 

7 3 = hyh* + hhz (3-33) 

h = hh-hz 2 (3.34) 

h = Ixh^ + hyhz (3.35) 

h = hh-hy 2 (3.36) 

7 7 = h — h ( 3 - 37 ) 

7g = I x — h (3.38) 

h = h~h (3.39) 
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In these equations, 7 Z , I y and I z are the moments of inertia about the body x-axis, 
y-axis and z-axis, respectively. Since the inertia matrix, 7 is symmetric, only three 
cross moments of inertia, 7 xy , I xz and I yz are required in addition to 7 X , l y and 
l z for the complete specification of the inertia matrix. However if the aircraft is 
symmetric about the x-z plane as is mostly the case, two components of the cross 
moments of inertia I xy and I yz can be assumed to be zero. This assumption will 
result in considerable simplification of the moment equations. 

Finally, the navigation equations are: 

ib U 

m =1^ v (3.40) 

it W 

The transformation matrix from the inertial to the body frame Tb/\ given by 
Equation (3.3). 

The six-degree-of-freedom model driven by the forces and moments is sum- 
marized in the block diagram given in Figure 3.4. It may be observed that the 
model requires twelve initial conditions. The only external inputs are the surface 
deflections and the throttle commands. 

The location of any light on the runway in the camera coordinates can be 
determined if the position and orientation of the aircraft and the coordinates of 
the runway lights with respect to the inertial frame are known. A mathematical 
model of the camera is needed to establish the relationship between the position 
of the runway lights in the camera frame, and their position in the image plane. 
Such a model will be described in the following section. 

3.4 Camera Model 

In order to avoid the complexities of having to deal with optical aberrations 
caused by lenses, it is customary [27] to represent the camera model by a pinhole 
lens together with an image plane located at the focus. The distance between the 
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Figure 3.4: Aircraft equations of motion. 
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lens and the focus is the focal length, /. The image of a point p in the scene 
is determined by a ray projected from the point through the lens center. The 
location where this ray intersects the image plane is where the image of the point 
is registered. Such a model is shown in Figure 3.5. The model in the figure results 



Figure 3.5: Pinhole camera model. 

in images that are inverted left to right and top to bottom. This is in contradiction 
with how 1 the human observer views the scene and how a television camera outputs 
the image. To avoid the inversion, a mathematically equivalent projection called 
the central projection [27] can be used. 

Central projection involves projecting a ray from a point to a frontal image 
plane such that the ray passes through the lens center. The geometry is shown in 
Figure 3.6. The camera axis system has its origin at c as previously shown in Figure 
3.2. The x-axis of the camera coordinate system is aligned with the optic axis. The 
z-axis points down and the y-axis completes the right-handed coordinate system. 
It can be observed that the central projection is a many-to-one mapping since 
all the object points along a projection ray are mapped to a single location in the 



74 



Frontal Image Plane 


Figure 3.6: Pinhole camera with frontal image plane. 

image plane. Thus if the location of an object point is known in the image plane, all 
that can be said about its three dimensional location is that it is located somewhere 
along the line passing through the image point and the lens center. On the other 
hand, if the location of the object point is known in the camera frame, its location 
in the image plane can be determined uniquely. This process is termed as direct 
perspective projection. These facts imply that there isn’t enough information in 
one image to recover three dimensional geometry. Two or more images obtained 
from different vantage points may be required to reconstruct the three dimensional 
scene. The process of recovering the three dimensional coordinates from one or 
more images is called inverse perspective projection. 

Real cameras capture the scene at discrete pixel locations indexed by rows 
and columns. Thus, every pixel is referenced by two coordinates: u and v with 
respect to a coordinate frame called the image frame with its origin, o, located at 
the top left hand corner of the image plane as shown in Figure 3.6. The u-axis 
is directed from left to right and the v-axis is directed from top to bottom of the 
image plane. Let, u c and v c be the coordinates of the camera center with respect 
to the image frame origin o. The image coordinates u p and v p for an object p 
can be obtained by constructing two sets of similar triangles from the geometric 
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relationships shown in Figure 3.6. Thus, 

u p — u c 

~T~~ 

Vp - Vc 

f 

These two equations describe the direct perspective projection process. The posi- 
tion components used here are defined in Equations (3.7) through (3.9). Note that 
the parameter, / is the focal length of the camera. 

Direct perspective projection equations for a pinhole camera model that take 
radial lens distortion, uncertainty scale factors and sampling into account are avail- 
able in the literature [98]. Such models are too complex for the purposes of this 
report. However, they may be useful for error analysis of cameras. 


yep 

X cv 

% Cp 


^cp 


(3.41) 

(3.42) 


3.5 Data Sources 

Pilots use airport lighting for obtaining alignment guidance and glide slope 
information during night approach and landing. The geometric information and 
color coding in the airport lighting layout is utilized by the human perceptual sys- 
tem for estimating position and orientation with respect to the runway. Since the 
human visual system sees a perspective image of the airport lighting, position and 
orientation estimation requires the pilot to correlate the scene with the aircraft 
position. In practice, this is accomplished by repeated landings at particular run- 
ways. Based on the parallels drawn between the human perceptual system and a 
conceptual machine vision system described in Section 1.6, it should be possible 
to derive the runway relative position and orientation information by comparing 
the airport images with the geometric model of the airport lighting layout. 

Before venturing into developing a system capable of generating the kind of 
information the pilot requires, an understanding of the standard airport lighting 
geometry is required. 
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3.5.1 Airport Lighting 

The purpose of airport lighting is to provide information about airport/runway 
identification, approach direction, alignment and attitude information for safe 
night landing. Standard airport lighting is composed of the approach and run- 
way lights [49]. The approach lights consist of centerline bars, sequenced flashers, 
threshold lights, cross bar lights, wing bar lights and the terminating bar lights. 
The runway lights include edge lights, centerline lights and touchdown zone lights. 
The approach as well as the runway lights are color coded and are located at fixed 
distances with respect to the runway threshold. 

In airports with several runways, the type of lighting used on the runways 
assists the pilot in determining if the aircraft is headed towards the desired airport 
and the correct runway. Sequenced flashers and color coding of the threshold bar 
indicate the approach direction to the pilot. Similarly the centerline and edge lights 
provide lateral and vertical alignment guidance. Additional information available 
in the airport lighting structure useful for a machine vision system is described in 
the following. 

3.5. 1.1 Standard Approach Lighting System 

Common configurations for approach lighting are the Calvert system and the 
standard configuration-A system [49]. The Calvert system is widely used in Europe 
and elsewhere around the world. The standard configuration-A approach lighting 
system is the national standard for civil and military use in the United States. 
Both systems are 3000 feet long. Figure 3.7 illustrates the standard configuration- 
A approach lighting system. 

In the standard configuration-A approach lighting system, the centerline bars 
are composed of five white lights separated by 40.5 inches. There is a sequenced 
flasher in front of each centerline bar. The distance between the centerline bars 
is 100 feet. The cross bar, located 1000 feet from the runway threshold, consists 
of eight white lights on each side of the centerline bar. These lights are separated 
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Figure 3.7: Standard configuration- A approach lighting system. 
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from each other by five feet. The threshold bar is located ten feet from the runway 
threshold and consists of green lights arranged five feet apart. The threshold 
bar runs along the entire width of the runway and extends 45 feet beyond the 
runway on each side. The wing bar is located at a distance of 100 feet from the 
threshold bar. The wing bar consists of five red lights placed symmetrically about 
the centerline. The inter-light separation is 40.5 inches. The terminating bar is 
located at a distance of 200 feet from the threshold bar, and consists of five red 
lights, 40.5 inches apart at the centerline and two sets of three red lights, five feet 
apart, placed symmetrically about the centerline. 

Approach lights are usually placed on pedestals of different heights. The 
specifications for mounting approach lights in the United States are available in 
Reference [29]. 

For operations in reduced visibility such as Category II or lower, the Interna- 
tional Civil Aviation Organization (ICAO) specifications are used for the lighting 
within 1000 feet of the runway threshold. The remaining 2000 feet of the lighting 
system is left as is. For the standard configuration-A system, this means that nine 
rows consisting of three red lights each are placed on either side of the centerline 
between the threshold bar and the the cross bar. Additionally, two rows with 
four white lights each are placed at 500 feet from the threshold, symmetrically 
about the centerline. Detailed layout of the Category II approach lighting system 
is described in Reference [49]. 

A medium approach lighting system (MALS) is often used at smaller airports 
for non-precision approaches. This system is 1400 feet long as opposed to 3000 
feet long standard configuration-A approach lighting system. Also, the threshold 
bar is not continuous. Only four lights are placed on each side of the threshold 
to indicate the indicate the approach direction. A MALS layout is also given in 
Reference [49]. 
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3. 5. 1.2 The Runway Lighting System 

The runway lighting system consists of edge lights, centerline lights and touch- 
down zone lights. A typical runway lighting layout is shown in Figure 3.8. Standards 



Figure 3.8: Runway lighting system. 

for design and installation of runway lighting systems are given in [33, 34]. 

The runway edge lights are high intensity white lights, except for the last 
2000 feet. The edge lights in the last 2000 feet are colored yellow to indicate a 
caution zone. The edge lights are located ten feet away from the pavement and 
the distance between the lights along-track is 200 feet. 

The centerline lights are 50 feet apart and run all the way to the end of the 




80 


runway. The first centerline light is 75 feet from the runway threshold. These 
lights are white, except for the last 3000 feet in the approach direction where they 
are color coded. The lights for the first 2000 feet are alternate red and white, while 
the last 1000 feet are all red. 

The touchdown zone lights start at 100 feet from the runway threshold and 
extend to 3000 feet in the direction of approach. The zone lights consist of three 
white lights, which are five feet apart and located at a distance of 30 feet about 
the centerline. The rows of zone lights are 100 feet apart from each other. 

3. 5. 1.3 Model of Airport Lighting 

In the two previous subsections, the layout geometry of the approach and 
runway lighting was described with respect to the threshold bar. Thus, by placing 
the origin of an inertial coordinate frame on the threshold bar with one axis aligned 
along the threshold bar and the other axis aligned with the runway centerline, the 
location of every light can be specified relative to the inertial coordinate system. 
These position coordinates form the airport lighting model. 

In order to construct the geometric model of the runway lighting, information 
available from a standard airport design text [49] was discussed in this chapter. 
The deviations from the standard layout for any airport in the United States are 
documented in Jeppesen Charts [54] and Instrument Flight Rules (IFR) Supple- 
ment [22], Detailed models of most airports can also be built using the information 
and survey maps available from city and county airport commissions. 

The airport lighting model provides one source of information for the machine 
vision systems. The images from an onboard camera forms the other source of 
information. 
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3.6 Landing Flight Trajectory and Image Se- 
quence Simulation 

The equations described in the earlier sections can be used for simulating 
the landing flight trajectories and the associated camera images of the airport 
lights. The purpose of such a simulation is to serve as a test bed for position and 
orientation determination algorithms developed in later chapters. 

The aircraft landing operation was described earlier in Section 3.1. A portion 
of this discussion was devoted to the desired glide path. The altitude, time and 
distance relationships for a approach speed of 220 feet/second along a three degree 
glide path were illustrated in Figure 3.1. 

These relationships can be used to obtain the conditions for simulation. These 
are: ( 1 ) the aircraft is initially at 400 feet altitude and 6633 feet downrange from the 
threshold, (2) the touchdown point is 1000 feet from the threshold, (3) the aircraft 
approach speed is 220 feet/second, (4) the aircraft sink rate is 11.5 feet/second and 
(5) the aircraft is perfectly aligned with the runway centerline. Thus, the descent 
path is given: 


x bc = —6633 + 220f 

(3.43) 

Vbc = 0 

(3.44) 

z bc = —400 T 11. 5f 

(3.45) 


where, t is the time from the initial position; x bc , y bc and z bc are the desired or 
commanded aircraft position components along the inertial x-axis, y-axis and the 
z-axis. 

Aircraft flight along the prescribed path can be simulated by using the air- 
craft aero-propulsive models along with the equations of motion discussed in a 
previous section. However, this requires aircraft specific aero-propulsive models 
and a suitable flight control system. These difficulties can be avoided by assuming 
that a suitable control system can be designed to closely track the trajectory. In 
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this case, the aircraft trajectory can be approximately simulated by using just the 
kinematic equations. Since the aircraft actual trajectory is expected to be close 
to the desired trajectory with minor deviations, the actual trajectory can be sim- 
ulated by driving the linear and angular acceleration components by white noise. 
Thus, the kinematic equations required for simulation are: 


Xb 


Vbx 

(3.46) 

Vb 

— 

V by 

(3.47) 


= 

Vbz 

(3.48) 

Vbx 

~ 

Vi b 

(3.49) 

Vby 

= 

Vib 

(3.50) 

Vbz 

= 

Vh 

(3.51) 

i' 


q sin 4> sec 6 + r cos (f) sec 6 

(3.52) 

6 

= 

q cos (f> — r sin <p 

(3.53) 

4> 

= 

p + q sin <f) tan 6 + r cos <j> tan 0 

(3.54) 

r 

= 

Vr 

(3.55) 

<7 

= 

Vi 

(3.56) 

P 

— 

Vp 

(3.57) 




(3.58) 


Here, U{, x , Vb y and Vb z are the components of the inertial velocity, V^; r/ Xb , r)y b and 
rji b are the white noise components driving the linear acceleration components; 
and rjr, r and 7/p are the white noise components driving the angular acceleration 
components. The nomenclature for other terms remain unchanged. 

The landing trajectory simulation is accomplished by integrating the system 
of Equations (3.46) through (3.58) with the initial conditions: Xb = —6633 feet, 
3/6 = 0 feet, Zb = —400 feet, Vb x = 220 feet/second, Vb y = 0 feet/second, Vb z = 11.5 
feet/second, ij> = 0 degrees, 6 = — 3 degrees, </> = 0 degrees, p = 0 degrees /second, 
<7 = 0 degrees/second, and r — 0 degrees/second. 
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The image generation process corresponding to the landing trajectory simu- 


lation is summarized in the block 


diagram in Figure 3.9. The camera position and 



WHITE NOISE 


Figure 3.9: Image generation process. 

orientation specified by the aircraft equations of motion, the camera model and a 
model of the airport lighting are used to generate an image of the airport lighting 
layout. 

Nighttime images of the airport are simulated using the lighting layout illus- 
trated in Figures 3.7 and 3.8. For image synthesis, the camera is assumed to be 
fixed to the aircraft, looking down along the glide slope. Since the camera axis is 
assumed to be colocated with the body axis, the look down angle is same as the 
pitch angle 0. The image is assumed to be digitized on a 512 x 512 pixel array, 
with the camera focal length being 600 pixels. This translates into a field-of-view 
of about 46 degrees. Image synthesis is achieved in two steps. First, the airport 
lights within the camera field-of-view are determined by using the known camera 
position, orientation and the field-of-view. Second, the lights within the field- 
of-view are projected onto the image plane by using the perspective projection 
equations described earlier. A simulated image constructed using this process is 
shown in Figure 3.10. This image corresponds to the camera located at an altitude 
of 95 feet and 812 feet downrange from runway threshold. 

The steps in the image generation process are summarized in Table 3.3. 
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Figure 3.10: Simulated airport image. 
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Table 3.3: Summary of the Image Generation Process 


1. Initialize the image matrix f(i, j) = 0; i = 1,2, . . . ,512; j — 1,2, . . . ,512. 

2. Use yaw, pitch and roll attitudes 6 and cj> to compute the elements of the 
transformation matrix r\ through rg using Equation (3.3). 

3. Compute the position of the airport lights x Cp , y Cp and z Cp using Equations 
(3.7), (3.8) and (3.9) with the actual aircraft position components, x b , y b 
and 26, such that p = 1,2, . . . , M where, M is the number of lights within 
the field-of-view r of the camera. 

4. Compute the location of each light in the image plane u p and v p using the 
perspective projection Equations (3.41) and (3.42). 

5. Quantize every u p and v p using [u p + 1/2J and [v p + 1/2J where, [ J 
the floor function. Following the definition of the floor function in Reference 
[38], |u p + 1/2J is the greatest integer smaller than or equal to u p + 1/2. 

6. Set the image matrix f{[u p + 1/2J, + 1/2J) = 256. 
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It should be noted that when the aircraft is at 400 feet altitude along three 
degree glide slope, both the approach lights and the runway lights are within 
the field-of-view. As the aircraft proceeds along the descent path, the approach 
lights begin moving out of field-of-view below 230 feet altitude. Any position and 
orientation determination algorithm must be able to adapt to this fact. 


3.7 Algorithm Development Considerations 

Algorithm development is dependent on the available data sources and math- 
ematical models. Two data sources were identified in previous sections. These are: 
(1) image of the airport lighting acquired by the camera and (2) the known airport 
lighting geometry. Earlier in this chapter, the airport light locations in the model 
were mathematically related to their respective locations in the image plane using 
a pinhole camera model and camera motion parameters. The image formation pro- 
cess using the mathematical relations was further discussed in Section 3.6. Earlier 
in Section 2.3.7, the difficulties of correlating the image features with the model 
features were examined. Clearly, these difficulties can be eliminated if the image 
features and model features are transferred into a common framework. There are 
two natural choices for performing such comparisons. The comparisons can be 
carried out in the inertial plane or in the image plane. Each of these choices result 
in different families of algorithms. 

Figure 3.11 illustrates the inertial frame-based family of methods for runway 
relative position and orientation. The image of the runway lighting acquired by 
the camera is transformed to the inertial plane using inverse perspective projec- 
tion. This requires an initial estimate of aircraft position and attitude. Equations 
(3.41) and (3.42) are used with Equations (3.7) through (3.9) to recover the inertial 
locations of the lights x p and y p by assuming a camera position, orientation and 
that all the airport lights are located on the z p = 0 plane. The difference between 
the features extracted from this layout and the features extracted from the known 
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Figure 3. IF Solution family I. 
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lighting layout form the inputs to the position and orientation algorithm. The 
algorithm iteratively improves the position and orientation estimates in order to 
achieve a better matching. The improved estimates are then used for inverse per- 
spective projection. This procedure recovers the camera position and orientation 
by driving the feature errors to zero. 

The procedure given in Figure 3.12 describes the second family of methods 
for runway relative position and orientation estimation by carrying out feature 
matching computations in the image plane. In this case, the camera model is 
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Figure 3.12: Solution family II. 

used together with the airport lighting model to predict the image of the lighting 
arrangement. This prediction is based on an assumed camera position and orienta- 
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tion. Next, the difference between the features extracted from the predicted image 
and those extracted from the actual image of the airport lighting as seen by the 
camera are fed into a position and orientation algorithm to refine the position and 
orientation estimates. These estimates are then used for updating the predicted 
image. Thus as in the first family of methods, the position and orientation states 
are recovered by driving the feature errors to zero. 

The two family of solution methods have their advantages and disadvantages. 
Since inverse perspective projection is used in the first family of methods, the as- 
sumption that all airport lights lie on the z p = 0 plane is necessary. Moreover, due 
to perspective projection, lights farther down the runway are bunched together in 
the image plane. Location of these lights with respect to the inertial frame cannot 
be accurately recovered using the inverse perspective projection. Additionally, this 
family of methods require active model adaptation to remove the lights outside the 
field-of-view of the model as the aircraft moves because portions of the airport 
along the descent path. The main advantage is that since the structure of the 
predicted and model lighting is well defined in the inertial frame, the comparisons 
are straight forward. 

The second family of methods use direct perspective projection to synthesize 
the predicted image. As a result, the assumption that all airport lights lie on 
the ground plane z p = 0 plane is not needed. The model adaptation process is 
automatic because only the lights that are within the field-of-view of the pinhole 
camera model are used for synthesis of the predicted image. The main disadvantage 
of this family of methods is that the structure of lighting in the predicted and the 
camera images is difficult to identify due to perspective distortion. 

Note that in both the procedures outlined in the foregoing, a single image of 
the airport lighting is used as a part of an iterative scheme to recover the camera 
position and orientation coordinates. Algorithms based on each of these families 
of solution methods will be discussed in the ensuing chapters. 
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3.8 Summary 

The foundations for the development of vision-based algorithms for position 
and orientation determination were laid in this chapter. First, the kinematics 
of flight along the glide slope were described. Federal Aviation Administration 
defined landing categories together with the associated decision heights, and the 
performance specification for precision approach and landing for these categories 
were then discussed. Since the glide slope is defined with respect to the runway 
centerline and the touchdown point on the runway, and the airport lighting ge- 
ometry was also defined with respect to the runway threshold. The origin of the 
inertial coordinate frame was assumed to be located at the threshold with the x- 
axis pointing along the centerline in the approach direction, the y-axis along the 
threshold and the z-axis pointing down. With this choice of the inertial coordinate 
system, and definitions of the body and camera coordinate systems, the locations 
of lights in the lighting model were related to their location with respect to the 
camera coordinate system. This relation was established in terms of the aircraft 
location and orientation with respect to the inertial coordinate system, and the 
camera location with respect to the aircraft fixed body coordinate system. Since, 
both the aircraft position and the orientation are a consequence of the transla- 
tional and rotational motion of the aircraft subjected to propulsive, gravitational 
and aerodynamic forces, the equations of motion describing the dynamics and the 
kinematics of the aircraft were discussed. A pinhole camera model was then de- 
scribed for relating the camera relative coordinates of the model lights to their 
image coordinates. Using the models and equations, procedures for simulating the 
landing flight and images along the landing path were described. Finally, the fact 
that direct and inverse transformations between the location of the lights in the 
inertial frame and the image plane can be computed using the equations described 
in the chapter resulted in two possible solution approaches for determination of 
runway relative aircraft position and orientation. 



Chapter 4 


Parameter Optimization Based 
Position Determination Methods 


The runway position determination techniques discussed in this chapter are 
based on the first solution family illustrated in Figure 3.11. The runway relative 
orientations, 6 and are assumed to be known in all the algorithms presented 
in this chapter. 

Let the image coordinates of a light p in the image be given by u p and v p . The 
relation between these coordinates and the camera relative position components is 
given by the direct perspective projection equations, described earlier in Equations 
(3.41) and (3.42). For notational simplification let, 


u P 

Vp 



v p - v c 


f 


(4.1) 

(4.2) 


where, u c and v c are the coordinates of the image center with respect to the image 


frame and / is the focal length of the camera. Substituting for u p and v p in terms 
of light and aircraft position vector components results in the following relations: 


r*{x E - Xb) + r 5 (y p - yk ) + r 6 (^ p - 2j,) + ky 

r x (x p - x b ) + r 2 {y p - y b ) + r 3 (z p - z b ) + k x 

rjjxp - Xh ) + r 8 (y p - y b ) + r 9 (^ p - z b ) + k z 

ri(x p x b ) -(- r 2 (jjp y b ) rz(zp z b ) -|- k x 
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Since the camera location with respect to aircraft center of gravity k x , k y and 
k z are known quantities, they can be set to zero without loss of generality. This 
simplification implies that the camera frame is colocated with the body frame. The 
inverse perspective projection equations can be obtained from Equations (4.3) and 
(4.4) as: 

z p = x b + ^^—^^ (z p -z b ) (4.5) 

*^lp^5p $2p$4p 

y, = y t + - Zb ) ( 4 . 6 ) 

^2p*^4p ^Ip^Sp 

The quantities s lp through se p depend on the image coordinates U v and V p , and 
the known elements of the transformation matrix rq through r 9 . They are defined 
by the equations: 


Sip - U p r\ - r 4 

(4.7) 

S 2 P = U p r 2 - r 5 

(4.8) 

■S3 p = Up r 3 — r 6 

(4.9) 

S4p = V P ri - r 7 

(4.10) 

$5p — Vpf2 ~ T 8 

(4.11) 

s 6p = Vp r 3 ~~ r 9 

(4.12) 


Examination of the inverse perspective projection equations given by Equa- 
tions (4.5) and (4.6), reveals that the x p and y p position components of all airport 
lights are shifted by the aircraft position components xi , and y b , and scaled by 
the aircraft altitude, — z b . Given the aircraft position components and the verti- 
cal coordinate z p of each light, Equations (4.5) and (4.6) uniquely determine the 
individual light horizontal position components x p and y p . In order to make the 
problem solvable, z p can be set to zero or a constant for all airport lights. This 
assumption is reasonable specially when the aircraft is at higher altitudes. Since 
the position of the aircraft relative to the plane containing the airport lights is of 
interest, define the aircraft altitude above the plane of the runway by h. Thus, 
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with z p = 0 assumption: 

x p = x 6 + g2pg6p ~^ 5p fe (4.13) 

p$5p $2p$4p 

V, = y„ + (4.14) 

^2p^4p ^lp^5p 

These two equations show that the projected positions of the lights on the plane 
of the runway are dependent on three parameters, z/& and h. 

Since the position of each airport light is known from the airport lighting 
layout, the position components x^, yb and h that results in a match between the 
known position and the projected position of each light is the desired solution. 
This problem can be conveniently posed as parameter optimization problem for 
determining aircraft position components that minimize the matching error. Early 
versions of two algorithms based on this idea are described in Reference [15]. Re- 
vised version of these algorithms are described in the following sections. The first 
algorithm assumes that the aircraft altitude is available from an onboard altimeter. 
The second algorithm does not make this assumption. 


4.1 Algorithm I 

This algorithm assumes that in addition to the aircraft attitudes, t/j, 6 and 
<^>, the aircraft altitude is available from an onboard altimeter. Thus, for every 
light detected in the image Si p through s 6 p can be computed using Equations 
(4.7) through (4.12). These can then be used in Equations (4.13) and (4.14) for 
computing the relative position components of each light, detected in the image, 
as follows: 

x v -x b = S ^flL h (4.15) 

-Slp-S5p — $2p-S4p 

^lp^Sp ^3p^4p 1 f . 1 n \ 

2/p — Vb = — — — -h (4.16) 

s 2pS4p — -SipSSp 

In the ideal case, the projected airport lighting layout would be bounded by a 
rectangle of the same dimensions as the rectangle bounding the model lighting 



94 


layout as shown in Figure 4.1. 



Figure 4.1: Envelope matching in the inertial frame. 

Using Equations (4.15) and (4.16), the minimum and the maximum relative 
coordinates of the projected lighting can be found to be x Pmin , x Pmax i J/p mm an< ^ 
y V max ' The max i ma and minima define the rectangle on the horizontal plane that 
encloses the projected lighting layout as shown in Figure 4.1. Since the coordinates 
of every airport light are known with respect to the inertial coordinate system, via 
the airport lighting model, the minimum and maximum coordinates of the lighting 
model can be found to be ar, min , x lmar , and yi max ■ These coordinates define 

the rectangle that encloses the model lighting shown in Figure 4.1. 

Consider the coordinates of the upper left corners, A’ and A of the enclosing 
rectangles in Figure 4.1. These are, {x Pmax , y Pmin ) with respect to the aircraft and 
(x tmax , yi m j n ) with respect to the inertial coordinate system. Using the relative 
geometry shown in Figure 4.1, it may be observed that: 

Xb = x, maI -x Pmax (4.17) 
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Vb Vimin yPmin 


(4.18) 


Similarly, using the coordinates of the lower right corners, B’ and B, 

%b ^tmtn ^Pmin (4.19) 

Vb = y^max-VPmax ^ 4 - 20 ^ 


The average of these two calculations can be taken as the aircraft position estimate. 
Combining Equations (4.17) with (4.19) and (4.18) with (4.20) yields: 


Xb = 


Vb 


( x imin X Pmin ) ( X *rnax X Pmax ) 

2 

(Z/imin ypmin ) ^TTiai ypmax ) 


(4.21) 

(4.22) 


So far only the ideal case has been examined. In the real situation, perspective 
projection causes the lights at the far end of the runway to be bunched together 
in the image plane making it difficult to determine their position. Since the lights 
closer to the camera are well separated in the image plane, it is more reasonable 
to use a weighted average of the x b position estimates. Thus an improved position 
estimate for Xb is: 

__ w l{ x tmin ~ X Pmin) max ~ X Pmax^ (4 23 ) 

Xb ~ (m + ^2) 

XD] and w 2 are the weighting factors. No such rationale can be applied along the 
y-axis. Hence, Equation (4.22) can be used directly. 

The position estimation algorithm is summarized in Table 4.1. 

The attractive features of this algorithm are: (1) envelope matching does not 
require any explicit identification of individual lights, (2) no iterative computations 
are required, and (3) only a single image is required. However, the two significant 
difficulties with this algorithm are that accurate knowledge of altitude is required 
and the elevation of the airport lights, — z p , is not included in the computations. 

Algorithm I does not take advantage of the fact that numerous images are 
available along the descent path. Since these images are related to each other by 


Table 4.1: Summary of Algorithm I 


1. For every light in the scene compute U p and V p using Equations (4.1) and 
(4.2) with p = 1,2, . . . , M where, M is the number of lights detected in the 
image. 

2. Us e ip, 0 and (f> to compute r\ through r 9 using Equation (3.3). 

3. Compute Si p through Se p using Equations (4.7) through (4.12). 

4. Compute x p — x;, and y p - Vb, P = 1,2, . . . , A/, using Equations (4.15) and 
(4.16) along with the known altitude. 

5. Compute mai(i p }, mm{x p }, max{y p } and min{y p }. 

6. Compute max{x,}, mm{x,}, max{yi] and i = 1 , 2 ,..., AT, using 

the coordinates of N airport lights within the model. 

7. Compute the inertial position components, Xf, and yt, using Equations (4.23) 
and (4.22). 
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aircraft motion, equations of motion can be used with this algorithm for obtaining 
improved aircraft position estimates. Kalman filter [3] is a natural choice for gener- 
ating improved position estimates by weighting the position propagated using the 
equations of motion and the position estimates provided by the algorithm in Table 

4.1 in a statistically optimal manner. For this purpose, a six-state Kalman filter 
with three position components and three velocity components is used together 
with Algorithm I for obtaining aircraft position and velocity estimates along the 
descent path. The details of this Kalman filter are described in Appendix A. The 
matrices required for implementing the Kalman filter are described in Appendix B. 
The two position components generated by Algorithm I are used as measurements 
for the Kalman filter. The known altitude is used as the third measurement. With 
these measurements, the Kalman filter provides the improved position and velocity 
estimates. 

Most genera] aviation aircraft use barometric altimeter which has a limited ac- 
curacy because the measurement depends on ambient temperature. Although the 
measurement accuracy is sufficient for maintaining the required vertical separation 
for Air Traffic Control, it is inadequate for operations close to the ground. If alti- 
tude could be computed reliably, it could be used for augmenting the barometric 
altimeter reading. 

An algorithm that does not require altitude measurements is described in Sec- 
tion 4.2. However, the lack of altitude information results in an iterative algorithm 
because three position components are estimated from two inverse perspective pro- 
jection equations. 

4.1.1 Results Using Algorithm I 

Results of two cases obtained using Algorithm I are described in this section. 
The first case is obtained using w x = 1 and w 2 = 0. The second case is obtained 
using the weighting factors uq = 1 and w 2 = 1 in Equation (4.23). In both cases, 
the landing scenario discussed in Section 3.6 is used. Aircraft landing flight trajec- 
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tory and the images along the glide slope are also simulated using the procedures 
described in Section 3.6. For initializing the Kalman filter, errors of 1000 feet in 
the along track position x and 100 feet in the cross-track position y i are assumed. 
The inertial velocity components were all initialized to zero. 

4. 1.1.1 Algorithm I with wi = 1 and iv 2 = 0 

The error residuals of the runway relative position and velocity components 
are described in this section. The error residual is defined as the difference between 
the value estimated by the Kalman filter and the true value. 

The along-track position error residual is shown in Figure 4.2. It may be seen 



Figure 4.2: Along-track position error using Algorithm I with tuj = 1 and w 2 — 0. 
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from the figure that the along-track position estimate converges to within ±100 
feet in less than one second. 

The cross-track position error residual presented in Figure 4.3 shows that the 



Figure 4.3: Cross-track position error using Algorithm I with rui = 1 and u >2 = 0. 

cross-track position estimate converges to within ±5 feet in less than one second. 

Figures 4.4 and 4.5 show the along-track velocity (t»6 x ) and cross-track velocity 
(ut y ) error residuals. It may be observed from the figures that the along-track 
velocity estimate settles to within ±10 feet/second in five seconds and the cross- 
track velocity settles to within ±5 feet/second in less than one second. 

The position error residuals for a range of altitudes corresponding to the FAA 
landing categories are summarized in Table 4.2. The cross-track position ( j/j>) 
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error residual is in feet. By comparing Table 4.2 to Table 3.2, it may be seen that 
Table 4.2: Algorithm I with = 1 and w 2 — 0 Results 


Category 

V b 

CAT I 

±0.54 

CAT II 

±0.22 

CAT Ilia 

±1.49 

CAT Illb k c 

±0.85 


Algorithm I along with the Kalman filter generates position estimates which meet 
the navigation accuracy requirements for all the three categories. 

4. 1.1. 2 Algorithm I with w\ = 1 and w 2 — 1 

In order to assess the benefit of matching both the near and the far ends of the 
envelopes shown in Figure 4.1, the error residuals of the runway relative position 
and velocity components obtained using Algorithm I with w\ = 1 and w 2 = 1 are 
examined in this section. 

Figure 4.6 shows the along-track position error residual. It may be seen from 
the figure that the along-track position estimation error continues to grow as a 
function of time. This is due to the inability of correctly recovering the position 
components of distant lights using the inverse perspective projection equations. 
The dimension of the predicted envelope along the viewing direction is therefore 
shorter than the model envelope. In an attempt to match the shorter envelope 
to the longer model envelope, the aircraft position estimate is erroneously esti- 
mated. These erroneous position estimates when provided as measurements to 
the Kalman filter results in incorrect position estimates. From Figure 3.11 it may 
be seen that since the position estimates are used for inverse perspective projec- 
tion, grossly incorrect estimates of the aircraft position components would lead to 
a gross mismatch between the predicted lighting and model layouts to an extent 
that subsequent correct recovery of the aircraft position components may not be 
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possible. 

The along-track position result described in Figure 4.6 under-predicts the 
true position which could cause the pilot to overshoot the touchdown point. A far 
more dangerous situation would arise if the along-track position estimates over- 
predicted the along-track position because that could cause the pilot to land short 
of the runway. 

The cross-track position error residual is shown in Figure 4.7. This figure 
shows that the cross-track position estimate converges to within ±5 feet in less 
than one second. 



Figure 4.7: Cross-track position error using Algorithm I with = 1 and w 2 = 1. 


The along-track velocity and cross-track velocity error residuals are illustrated 
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in Figures 4.8 and 4.9. Figure 4.8 shows that the along-track velocity estimation 
error never settles down. The cross-track velocity estimates on the other hand 
settle to within ±5 feet /second in less than one second. 



Figure 4.8: Along-track velocity error using Algorithm I with Wi - 1 and w 2 — 1- 

The cross-track position error residuals for FAA landing categories are summa- 
rized in Table 4.3. It can be observed that the navigation accuracy requirements 
for all the three categories are met. Note navigation accuracies for the landing 
categories are not specified for the along-track position. 

Comparing the results of Subsections 4. 1.1.1 and 4. 1.1. 2 suggests the weights 
Wi = 1 and w 2 = 0 yields better position estimates. 


Cross-track Vel. Error Res. (ft/s 
— 1 o -5 0 5 10 
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Figure 4.9: Cross-track velocity error using Algorithm I with u>i — 1 and — 1. 


Table 4.3: Algorithm I with tui = 1 and u> 2 = 1 Results 


Category 

Lateral (yt) 

CAT I 

±0.97 

CAT II 

±0.53 

CAT Ilia 

±3.73 

CAT Illb k c 

±0.44 
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4.2 Algorithm II 


The main difference between this algorithm and the previous one is that the 
present algorithm estimates all the three position components of the aircraft. 

As in the previous algorithm, the minimum and maximum values timin') 
xii and Vi „ can be determined using the known x, and t/, coordinates of the 
lights in the model. However, the lack of knowledge about the aircraft altitude 
makes the process of projecting the lights detected in the image onto the horizontal 
plane somewhat more complicated. The projection requires the assumed values of 
aircraft position components, x b , yb and h. The quantities through Se p in these 
equations can be evaluated in exactly the same manner as in Algorithm I. The 
position coordinates x v and y p of every light can be found using the guessed initial 
position. These coordinates can then be used to find the values x Pmin , x Vmax , yp m i n 
and y Pmax - Since the rectangle formed by the maximum and minimum values 
are required to enclose the projected and model lighting layouts, the position 
determination problem can be viewed as a parameter optimization problem for 
matching the projected rectangle with the rectangle enclosing the airport light 
database. 

A quadratic cost function for measuring the matching error can be constructed 
in terms of the coordinates of the upper left and lower right corners of the enclos- 
ing rectangles, (^p max i 2/p m i n )» (^'Pmin 5 yPmax ) ^ 

follows: 


J = 


+ 


w l { x imi n ^Pmin) x Pmax ) 

) 2 




{w\ + w\) 

(y*min ~ VPminY + (V imax yPmax > 


sft 


(4.24) 


As in Algorithm I, the weights 0 < W\ <1 and 0 < w? < 1 can be used to 
establish the error contribution of lights near the camera and those that are far 
away. The objective of the optimization algorithm is to determine aircraft position 
components x b , y b and h that minimize the performance index J. Any one of the 
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several unconstrained optimization methods can be used for this purpose. The 
Newton-Raphson method is chosen for the present research. 

The Newton-Raphson method [74] is an iterative method for finding the zeros 
of a function. In one dimension, the Newton-Raphson technique extrapolates the 
derivative at the current location in order to find an improved estimate of the 
zero. The method has its basis in the Taylor Series expansion of a function in the 
neighborhood of a point. It is known to converge quadratically when the function 
is smooth and convex. The Newton-Raphson technique can be readily extended to 
multiple dimensions. The technique is suitable for solution of nonlinear systems of 
equations. The technique can be adapted for finding the extremum of a function 
by driving the gradient vector of the function to zero. The difference between 
searching for a zero of a function and the zero of a derivative is that, in the first 
case the Jacobian matrix is used while in the second case the Hessian matrix is 
used. 

Thus, the Newton-Raphson formula for position determination is: 


£ b 


Xb 
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Vb 
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y b 

- H - 1 

dJ 
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QdIQj 
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(4.25) 


where, the Hessian matrix H is given by: 
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d 2 J 
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dhdxb 


d 2 J 
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dhdyb 

d 2 J 

d 7 J 

d 2 J 

dxbdh. 

dybdh 

dh 2 


(4.26) 


The subscript n + 1 denotes the improved estimate of the aircraft inertial position. 
Starting from an initial guess, the position vector can be iteratively computed 
using Equation (4.25) until the change in the cost function is smaller than a preset 
tolerance. 

Equation (4.25) requires the computation of the first and second partial deriva- 
tives of the cost function. Since the cost function depends on ( 1 Pmax’ VPmin) an< ^ 
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v min' V? max) coordinates, their partial derivatives are required for computing the 
partial derivatives of the cost function with respect to aircraft position coordinates. 
These can be derived from Equations (4.13) and (4.14) as follows: 
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(4.27) 

(4.28) 

(4.29) 

(4.30) 

(4.31) 

(4.32) 

(4.33) 

(4.34) 

(4.35) 

(4.36) 

(4.37) 

(4.38) 


Note that these partial derivatives are evaluated at the instantaneous values of 
Xr, , x„ . , u„ and y v . Partial derivatives of the cost function with re- 
spect to aircraft position vector components can be computed analytically by using 
Equation (4.24) and Equations (4.27) through (4.38) as: 


dJ_ 
dx b 

dJ_ 

dy b 


|tUl(Xj m i n Xp min ) + W2(Xi ma x I Pmax)} 


sj{w\ + u|) 
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(4.39) 

(4.40) 
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The elements of the Hessian matrix are given by: 
d 2 J 2(w\ + w 2 ) 
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(4.41) 

(4.42) 

(4.43) 

(4.44) 

(4.45) 

(4.46) 

(4.47) 

(4.48) 

(4.49) 


(4.50) 


Since the gradient vector and Hessian computations are analytic, the Newton- 
Raphson iterations can be carried out at a high computational rate. The essential 
steps involved in Algorithm II are summarized in Table 4.4. 

The three aircraft position components estimated using Algorithm II are used 
as the measurements of the six-state position/velocity Kalman filter. As in Algo- 
rithm I, the Kalman filter integrates the information derived from multiple images 
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Table 4.4: Summary of Algorithm II 


1. Compute max{y,} and min{yi } using the coordinates 

of N airport lights within the model. Here, 1 < i < N. 

2. Compute U p and V p V p, such that 1 < p < M where, M is the number of 
lights detected in the image, using Equations (4.1) and (4.2). 

3. Use tp, 0 and cf> to compute rq through r 9 using Equation (3.3). 

4. Compute s ip through s 6p using Equations (4.7) through (4.12). 

5. Assume Xb, yt and h. 

6. Using Equations (4.13) and (4.14), compute x p and y p V p, such that 1 < 

p < M . 

7. Compute max{x p ), min{x p }, max{y p } and min{y p }. 

8. Compute cost J 2 using Equation (4.24). 

9. If this is the first time, skip to step 11; else, continue. 

10. Is | J 2 - Ji |< 8 ? If yes, stop; else, continue. The parameter 8 is the 
stopping tolerance. 

11. Set Ji = J 2 - 

12. Compute the partial derivatives using Equations (4.39) through (4.50). 

13. Using the Newton-Raphson formula, Equation (4.25), compute the inertial 
position components, Xb and t/j, and h. 

14. Return to step 6. 
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processed using Algorithm II. The six— state Kalman filter used here is described 
in Appendix A. The matrices required for this Kalman filter are available in Ap- 
pendix B. 

4.2.1 Results Using Algorithm II 

Results of the two cases obtained using Algorithm II are described in this 
section. The first case is obtained using w\ = \ and w 2 = 0 while the second case 
is obtained using uq = 1 and w 2 = 1 in Equation (4.24). As in Algorithm I, the 
landing scenario discussed in Section 3.6 is used for both the cases. Errors of 1000 
feet in the along-track position x b , 100 feet in the cross-track position y b and 100 
feet in the altitude h are assumed for initializing the Kalman filter. The inertial 
velocity components are all initialized to zero. 

4. 2. 1.1 Algorithm II with ttq = 1 and w 2 = 0 

The along-track position error residual is shown in Figure 4.10. It may be 
seen from the figure that the along-track position estimate converges to within 
±100 feet in less than one second. 

The cross-track position error residual presented in Figure 4.11 shows that 
the cross-track position estimates converge to within ±5 feet in less than three 
seconds. 

The altitude error residual is shown in Figure 4.12. It may be observed from 
the figure that the altitude error converges to ±5 feet within two seconds and 
stays within these bounds up to 30 seconds. Beyond that point the altitude error 
increases. The error increase can be attributed to the reduction in the number of 
lights within the field— of— view. The altitude of the aircraft at 30 seconds is 55 feet. 

Figures 4.13, 4.14 and 4.15 show the along-track velocity (u bl ), the cross-track 
velocity (u fcv ) and the sink rate (-v bz ) error residuals. It may be observed from 
the figures that the along-track velocity estimate settles to within ±10 feet/second 
in less than five seconds, the cross-track velocity settles to within ±5 feet/second 
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in less than one second and the sink rate settles to within ±5 feet/second in less 
than two seconds. These results are comparable to those obtained using Algorithm 
I for the same set of weights. 

Table 4.5 lists the position error residuals. In Table 4.5, the position com- 
ponents are in feet. Comparing Table 4.5 with Table 3.2, it may be seen that 

Table 4.5: Algorithm II with uq = 1 and w 2 = 0 Results 


Category 

Lateral ( yb ) 

Vertical ( h ) 

CAT I 

±0.70 


CAT II 

±0.22 


CAT Ilia 

±2.30 


CAT Illb k c 

±0.55 

±10.11 


Algorithm II along with the Kalman filter generates position estimates which meet 
the navigation accuracy requirements for Categories I and II. 

4. 2. 1.2 Algorithm II with w\ = 1 and w 2 = 1 

The error residuals of the runway relative position and velocity components 
obtained using Algorithm II with Wi = 1 and w 2 = 1 are described in this section. 

Figure 4.16 shows that the along-track position estimates converge to within 
±100 feet in about six seconds. 

The cross-track position error residual shown in Figure 4.17 shows that the 
cross-track position estimates converge to within ±5 feet in less than one second. 

The altitude error residual given in Figure 4.18 shows that the altitude esti- 
mates converge to within ±20 feet in less than six seconds. 

By comparing Figure 4.16 with Figure 4.6 it may be noted that the along- 
track position is significantly improved by Algorithm II at the cost of increased 
altitude estimation error shown in Figure 4.18. 

The along-track velocity, cross-track velocity and sink rate error residuals are 
illustrated in Figures 4.19, 4.20 and 4.21. Figure 4.19 shows that the along- 
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track velocity settles down to within ±10 feet/second in less than five seconds. 
The cross-track velocity estimates settle to within ±5 feet /second in less than one 
second. The sink rate estimates settle to within ±5 feet/second in less than five 
seconds. 

The position error residuals for ranges of altitudes are summarized in Table 
4.6. The nomenclature for this table is same as that used in Table 4.5. Table 4.6 

Table 4.6: Algorithm II with w\ = 1 and w<i = 1 Results 


Category 

Lateral (yt,) 

Vertical ( h ) 

CAT I 

±1.21 

±16.11 

CAT II 

±0.19 

±15.50 

CAT Ilia 

±3.52 

±10.72 

CAT Illb & c 

±0.62 

±9.01 


shows that in this case, the navigation accuracy requirements listed in Table 3.2 
are not met for any of the categories because of the reduced altitude estimation 
accuracy. 

As in Algorithm I, the results of these two cases indicate that Algorithm II 
should be used with = 1 and tc 2 — 0 weight combination. 


4.3 Summary 

Two methods for aircraft position estimation were described in this chapter. 
Both are based on envelope matching in the inertial frame. The first method re- 
quired the altitude in addition to the yaw, pitch and roll orientation angles. The 
second method did not require knowledge of aircraft altitude. Both the algorithms 
were based on the inverse perspective projection equations which relate the image 
coordinates of the airport lights to their inertial coordinates via the aircraft iner- 
tial position components. The problem of position determination was posed as a 
problem of parameter optimization with the aircraft inertial position components 
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as the parameters. In Algorithm I, the two position components were determined 
by matching the envelope of the predicted lighting geometry with the envelope of 
the stored lighting geometry model. In Algorithm II, a quadratic cost function 
was used for minimizing the envelope matching error using a Newton-Raphson 
method. The position components computed by these algorithms were used as 
measurements for a six-state Kalman filter for estimating the position and veloc- 
ity components along the glide slope. Results obtained for the four cases were 
discussed. It was shown that both Algorithm I and Algorithm II generate more 
accurate estimates with w-i = 1 and w 2 = 0 weight combination. This implies that 
these image based algorithms should place a higher weight on nearer lights, than 
on lights farther away along-track. 



Chapter 5 


Feature Correspondence Based 
Aircraft Position Estimation 
Methods 


Two position determination methods were presented in the previous chapter. 
The central idea there was to match the observed lighting layout to the model of 
the airport lighting layout. One of the difficulties with the algorithms described in 
Chapter 4 is that the cost function is based only on the coordinates of two corners 
of the envelope of the airport lights. Thus, they do not exploit all the information 
available in the image. The objective of the methods presented in this chapter is 
to formulate methods that use most of the information available in the image. 

Reason for using information from many lights as opposed to a few is obvious. 
An algorithm using information from multiple lights can be expected to be less 
sensitive to camera induced errors. Additionally, the estimation accuracy will 
remain unaffected even if some of the lights have failed. Algorithms presented in 
this chapter do not require any iterative computations. Significant benefits of such 
direct schemes are that they are computationally efficient and robust. 

The main idea in this chapter is to synthesize measures that capture the 
structure of the airport lighting layout in terms of scalar functions. Such measures 
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can be used as features for establishing the correspondence between the observed 
lighting layout and the model lighting layout. The important structures in the 
standard airport lighting layout were presented in Figures 3.7 and 3.8. Of these, 
the prominent ones are centerline, left and right edge lights and the threshold bar. 
It is possible to model these structures as straight lines using the Hough transform 
method discussed in Subsection 2.3.1. However, the Hough transform requires 
careful threshold selection to prevent phantom lines from being detected. Another 
approach for discovering the structure is to simply arrange the coordinates of the 
airport lights in a non-decreasing order. This arrangement causes the lights to be 
re-indexed so that their inertial coordinates are in a non-decreasing order. Any 
efficient method can be used for arranging. The Quicksort method [56] has been 
used in this research. For example, Figure 5.1 shows the y-coordinates of the 
airport lights arranged in a non-decreasing order. 

The basis for arranging the y-coordinates in a non- decreasing order lies in the 
fact that the left edge lights have a y-coordinate value of -100 feet, the centerline 
lights have a y-coordinate of zero feet and the right edge lights have a y-coordinate 
of 100 feet as may be seen in Figures 3.7 and 3.8. The graph shown in Figure 5.1 
is obtained by plotting the inertial y-coordinates of the re-arranged lights against 
their scaled indices. The scaled index for a light is obtained by normalizing its 
re-ordered index by the total number of lights in the model and multiplying the 
result by 100. Note that the scaled index is no longer an integer, but a rational 
number. The multiplication factor of 100 is chosen for convenience. This choice 
of multiplication factor results in the left edge lights having scaled indices with 
values less than 15 as may be seen in Figure 5.1. The range of scaled indices for 
the centerline lights and the right edges lights are also marked in Figure 5.1. The 
benefit of scaling is that the observed lighting layout is made comparable to the 
model lighting layout. Note that usually fewer lights are observed in the image 
compared to those available in the airport lighting model. 

The inertial along-track x-coordinates of the model lights arranged in a non- 
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decreasing order are shown in Figure 5.2. From the standard approach and runway 


* 



Figure 5.2: Along-track x-coordinates of the airport lights in a non-decreasing 
order. 

lighting layouts shown in Figures 3.7 and 3.8, it may be observed that the x- 
coordinates of the threshold bar lights have a common value of zero. Figure 5.2 
shows that threshold lights have scaled indices between 20 and 30. Other structures 
are difficult to identify in this graph because the x-coordinates of structures like 
centerline lights and edge lights are mixed together along the length of the runway. 
These structures appear prominently when their y-coordinates are arranged as 
shown in Figure 5.1. 

Since the observed and model lighting have the same range of scaled indices, 


132 


the corresponding structures are related by the range of indices. For example, the 
left edge lights in both the observed and model lighting layouts are expected to 
have scaled indices with values between 0 and 15. Thus, this range of scaled indices 
can be used for establishing the correspondence between features based on the left 
edge lights identified in both the layouts. 

Two position determination methods that exploit the ability to identify struc- 
tures in the graphs such as those in Figures 5.1 and 5.2 are described in this 
chapter. Both the algorithms belong to the first family of solution methods de- 
scribed in Figure 3.11. The first algorithm assumes that the altitude is known 
while the second algorithm does not make this assumption. Both algorithms as- 
sume that the aircraft attitudes, t/>, 0 and (f>, are known. Starting points for both 
the algorithms are the relationship between the aircraft position coordinates and 
the position of the lights given by the Equations (4.13) and (4.14). Thus, like the 
earlier algorithms, they are not designed to account for the varying elevation of the 
airport lights. The details of these algorithms are given in the following sections. 


5.1 Algorithm III 


Algorithm III assumes that in addition to the aircraft attitudes, tf>, 6 and 
4>, the aircraft altitude is available from an onboard altimeter. This assumption 
results in decoupling of Equations (4.13) and (4.14). Two such equations can be 
written for each light observed in the image. Adding the resulting expressions: 
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Here M is the number of lights. 


In the ideal case, the location of the lights in the 


airport lighting model x, and j/j should be same as that observed. That is: 
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E Vi = E Vv ( 5 - 4 ) 

l<><Af 1 <p<M 

Thus, 

E x, = Mx. + A ■£ I * 3|, * 5f ! < 5 ' 5 > 

1 <i<M l<p<M ( S lp S 5p s 4p$2p) 

E * = m» + '> E I (5-6) 

l<t<M l<p<M V 5 2 p 54p -5ipS5 p j 

In reality, the left hand sides are known exactly from the airport lighting 
model, while the right hand sides are not known exactly due to observation errors 
in the image. For example, perspective projection causes distant edge lights and 
centerline lights merge together to form lines in the image in Figure 2.2. Due to 
this and other effects, the exact location of lights in the image plane cannot be 

accurately ascertained. In order to make the left and right sides comparable, the 

left hand side is normalized with the number of N lights in the airport lighting 
model, while the right hand sides are normalized using the number of M distinct 
lights observed in the image. Thus, 
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In this form, the inertial aircraft position components and t/f, can be computed 
simply by subtracting the image based arithmetic means from the model based 
arithmetic means. 

Results obtained using Algorithm I and II in Chapter 4 suggested that only the 
nearer lights should be used for obtaining the along-track position of the aircraft. 
Therefore, only N\ nearer model lights and Mi nearer image lights should be used 
in Equation (5.7). Note that Ni < N and Mi < M. Using the nearer lights with 
scaled indices between zero and 40, Equation (5.7) can be re-written as: 
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It may be observed from Figure 5.2 that the approach lights and the threshold 
bar lights in the model have scaled indices between zero and 40. The cross-track 
component of the aircraft position can be obtained by using all the model and 
image lights in Equation (5.8). 

An additional fact that deserves consideration is that only a portion of the 
airport lighting is visible within the camera field-of-view as the aircraft follows 
its descent trajectory. Thus, the comparison between model lighting and actual 
lighting observed in the image may yield incorrect position estimates. In order to 
partially offset the errors caused by this, the position estimation algorithm employs 
the previously estimated aircraft position to determine the model lights within the 
field-of-view. This information is then used in the calculations using Equations 
(5.8) and (5.9). 

This does not cause any problems for initialization because the complete air- 
port lighting is visible to the camera at the beginning of the descent path even 
with substantial differences between the actual and assumed aircraft positions. 

The modified airport lighting model x and y coordinates arranged in a non- 
decreasing order are shown in Figures 5.3 and 5.4. These figures were generated 
using an assumed along-track position which was in error by 1000 feet compared 
to the true aircraft position. No errors were assumed in the cross-track position 
and altitude. It may be observed from the figures that the threshold bar lights, 
edge lights and the centerline lights are correctly represented in the model layout 
graphs. The figures also show the observed lighting x and y coordinates arranged in 
a non -decreasing order. Since the observed lighting is with respect to the aircraft, 
the true along-track and cross-track positions of the aircraft can be obtained by 
shifting the observed layout graphs to the model layout graphs. Figures 5.3 and 
5.4 show that the true along-track and cross-track positions are -6633 feet and 0 
feet with respect to the inertial coordinate system. 

Algorithm III is summarized in Table 5.1. 

Some of the salient features of Algorithm III are as follows. Firstly, since 
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Figure 5.3: x-coordinates of adapted model and image lights in a non-decreasing 
order. 
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Table 5.1: Summary of Algorithm III 


1. Arrange x, coordinates of all the N lights from the airport lighting geometry 
model in non-decreasing order. 

2. Compute the scaled indices k t = \00i/N] 1 < i < N. 

3. Select the ordered i, V i such that 0 < fc, < 40. Compute the mean of the 
N\ selected X; in Equation (5.9). 

4. Compute the mean using all y t , 1 < i < N, in Equation (5.8). 

5. For every light in the scene compute U p and V p using Equations (4.1) and 
(4.2) with p = 1,2, . . . , M where, M is the number of lights detected in the 
image. 

6. Use 0, 9 and <f> to compute rq through r 9 using Equation (3.3). 

7. Compute s lp through s$ p using Equations (4.7) through (4.12). 

8. Compute X p = (S2 p 56p — S 3 P S 5 P )/{ S 1 P S5 P ~~ ■ S 4 P S 2 P ) 

and y p = (si p s 6p - s 3p s ip )/(s 2p s 4p - si p s 5p ) Vp; 1 < p < M. 

9. Arrange x p coordinates of M image lights in non-decreasing order. 

10. Compute the scaled indices k p = lOOp/M; 1 < p < M. 

11. Select the ordered x p V p such that 0 < k p < 40. Compute the mean of the 
Mi selected x p in Equation (5.9). 

12. Compute the mean using all p p , 1 < p < M, in Equation (5.8). 

13. Use the known altitude h in Equations (5.9) and (5.8) and subtract the 
image-based means from the model-based means in these equations to com- 
pute the aircraft position components x*, and t/*,. 
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the method uses information from multiple airport lights, it is more robust when 
compared with the envelope matching methods discussed in the last chapter. Air- 
craft position is determined without resorting to the traditional correspondence 
methods such as the local correlation methods. 

Major limitations of the algorithm are that it is unable to account for the 
runway elevation and that knowledge of aircraft altitude is required. Due to the 
importance of accurate altitude information during landing, it is desirable to mod- 
ify this algorithm to provide the altitude information in addition to the lateral and 
longitudinal positions. This is the motivation for the next algorithm for aircraft 
position determination described in Section 5.2. 

As in the case of Algorithms I and II, the position measurements from individ- 
ual frames can be made consistent with the aircraft kinematics using a six-state 
Kalman filter described in Appendix A. 

5.1.1 Results Using Algorithm III 

The position and velocity estimates generated using the six-state Kalman 
filter driven by the measurements from Algorithm III are described in this section. 
These estimates were obtained using the aircraft landing flight trajectory and the 
images along the glide slope simulated using the procedures described in Section 
3.6. Errors of 1000 feet in the along-track position and 100 feet in the cross- 
track position are assumed for initializing the Kalman filter. The inertial velocity 
components are all initialized to zero. 

The along-track position error residual is shown in Figure 5.5. It may be seen 
from the figure that the along-track position estimate converges to within ±100 
feet in less than one second. 

The cross-track position error residual presented in Figure 5.6 shows that the 
cross-track position estimate converges to within ±5 feet in less than one second. 

The along-track velocity error residual is shown in Figure 5.7 and the cross- 
track velocity error residual is shown in Figure 5.8. It may be observed from the 
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figures that the along-track velocity estimate settles to within ±10 feet/second in 
less than six seconds and the cross-track velocity settles to within ±5 feet/second 
in less than one second. 

The position and velocity error residuals are summarized in Table 5.2. In 
this table the lateral component is in the units of feet. Comparing Table 5.2 


Table 5.2: Algorithm III Results 


Category 

Lateral (yt>) 

CAT I 

±0.61 

CAT II 


CAT Ilia 


CAT Illb & c 



to Table 3.2, it may be seen that Algorithm III along with the Kalman filter 
generates position estimates which meet the navigation accuracy requirements for 
all Categories. It may be seen by comparing Tables 4.2, 4.3 and 5.2 that Algorithm 
III results are of comparable accuracy as Algorithm I with wi = 1 and w 2 = 0. 
The results are better for Category III. 


5.2 Algorithm IV 

This algorithm is motivated by the desire to estimate the three position com- 
ponents of the aircraft relative to the inertial coordinate system located on the 
runway. In Chapter 4 the three position components were estimated using Al- 
gorithm II which was based on an iterative scheme. The attempt here is to use 
the structure discovered in Algorithm III to compute all three components of the 
aircraft position without the use of an iterative scheme. 

As discussed earlier, if any one of the aircraft inertial position components is 
known, the remaining position components can be recovered using Equations (5.8) 
and (5.9). However, in order to solve for three position components, at least one 
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additional equation relating the image quantities to the model quantities in terms 
of the aircraft position components is required. 

Several such equations can be formed by using the two basic Equations (5.7) 
and (5.8) on groups of lights that represent well defined structures in the ordered 
graphs. For example, an additional equation can be constructed from Equation 
(5.8) as follows. Scaled indices between zero and 15 of the ordered y-coordinates 
of the model lights and the observed lights can be used in Equation (5.8) to yield: 


W E Vi = Vb 

iV 2 l<i<JV 2 


+ 


h 

W 2 


E 

l<p<A / 2 


(•Slp^6p ^3p^4p ) 
(^’2p'®4p ^lp*^5p) 


( 5 . 10 ) 


Here, N 2 and M 2 are the number of lights that represent the left edge lights in the 
model lighting layout and in the observed lighting layout. 

A similar equation can also be written for the right edge lights by using the 
scaled indices between 85 and 100 of the ordered sets as: 


1 _ . h ifipfep ~ ■ S 3p-S4p) 

(N - N 3 + 1 ) N £ iN y '~ yh * M - M 3 + 1 M Jr<« - a,,.,,) 1 ’ 

where, N 3 and M 3 are the number of scaled indices between zero and 85 of the 
ordered sets. 

Re-writing Equations (5.9), (5.8), (5.10) and (5.11) as: 


Xb + CL\h 

= 6 X 

(5.12) 

Vb + a ih 


(5.13) 

Vb + a 3 h 

= ^3 

(5.14) 

Vb + ^4 h 

= b 4 

(5.15) 


it may be seen that Equation (5.12) and any one of the remaining three equations 
form a linearly independent set of two equations. Note, a x to a 4 are the image based 
averages and 6 X to b 4 are the model based averages. Equations (5.14) and (5.15) are 
linearly dependent if a 3 = a 4 . This can only happen if all the lights being consid- 
ered lie along a single projection ray from the lens center. Such a situation cannot 
arise when the camera is above the plane of the runway as discussed in Section 
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2. 3. 7.1. Therefore, Equations (5.14) and (5.15) are always linearly independent 
for the viewing geometry considered here. Equation (5.13) can be obtained by a 
linear combination of Equations (5.14) and (5.15) with a perfect camera. This is 
because the left and right edge lights are symmetrically located about the cen- 
terline lights. Errors in the imaging process can cause Equation (5.13) also to be 
linearly independent. 

Since Equations (5.14) and (5.15) are linearly independent, the aircraft alti- 
tude can be estimated by subtracting Equation (5.11) from Equation (5.10). Thus, 

^ _ N2 yj_ ~ (jv-N 3 +i) ^N 3 <i<N 2/« _ 10 , 

1 Y' ( 3 lp^6p ^3 p^4 p ) 1 Y' (^lp36p ^3p^4p) 

M2 2 ^ 1 <p<M 2 (S2p»4p — »lp*5p) M-M3+I 2^M3 <p<M (s 2p »4p— Slp^Sp) 

With the altitude so determined, Algorithm III can be used for determining the 
Xb and yb components of the aircraft inertial position. 

Alternatively, Equations (5.9), (5.8), (5.10) and (5.11) can be used for obtain- 
ing a least squares solution for the aircraft inertial position components. Following 
the notation of Equations (5.12) through (5.15), the resulting Least Squares solu- 
tion [55] is given by: 

[xf „yb,h] T = (A t A)~ l A T [b 1 ,b 2 ,b 3 ,b 4 ] T (5.17) 


where, 

1 0 ai 

0 1 a 2 

A = 

0 1 a 3 
0 1 a 4 


(5.18) 


It is important to point out that the proposed method estimates all the three 
position components without resorting to any iterative calculations. Thus, this 
approach can be expected to be faster and more robust than the Algorithm II 
described in Chapter 4. It may be noted that as in Algorithm III, the six-state 
Kalman filter driven by the outputs generated by Algorithm IV can be used for 
estimating the aircraft position and velocity components. 

Algorithm IV is summarized in Table 5.3. 
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Table 5.3: Summary of Algorithm IV 


1. Arrange and y, coordinates of all the N lights from the airport lighting 
geometry model in non-decreasing order. 

2. Compute the scaled indices ki = 100 i/N\ \ < i < N . 

3. Select the ordered x, V i such that 0 < fc, < 40. Compute the mean 6i for 
Equation (5.12) using the N\ selected x, in Equation (5.9). 

4. Compute 62 for Equation (5.13) using all y,, 1 < i < N, in Equation (5.8). 
Compute 63 and b 4 for Equations (5.14) and (5.15) by selecting the ordered 
y, such that 0 < &, < 15 and 85 < ki < 100 respectively and using them in 
Equations (5.10) and (5.11). 

5. Compute s \ p through se p using steps 5, 6 and 7 of Algorithm III in Table 
5.1. 

6. Compute x p = ( s 2p s 6p - s 3p s 5p )/(si p s 5p - s 4p s 2p ) 

and y p — (-Sip-S6p *^3p'^4p)/('^2p^4p *^ip^5p) Vp, 1 — P — A/. 

7. Arrange x p and y p coordinates of M image lights in non-decreasing order. 

8. Compute the scaled indices k p = lOOp/M; 1 < p < M . 

9. Select the ordered x p V p such that 0 < k p < 40. Compute the mean a j for 
Equation (5.12) using the Mj selected x p in Equation (5.9). 

10. Compute a 2 for Equation (5.13) using all y p , 1 < i < M, in Equation (5.8). 
Compute a 3 and a 4 for Equations (5.14) and (5.15) by selecting the ordered 
j/p such that 0 < k p < 15 and 85 < k p < 100 respectively and using them in 
Equations (5.10) and (5.11). 

11. Use Equation (5.17) to compute the aircraft position components x;,, ?/6 and 

h. 
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5.2.1 Results Using Algorithm IV 

The position and velocity estimates generated by the six-state Kalman filter 
driven by outputs of Algorithm IV are described in this section. The simulation 
scenario and the initial conditions were same as in Algorithm II in Chapter 4. 

The along-track position error residual is shown in Figure 5.9. It may be seen 



Figure 5.9: Along-track position error using Algorithm IV. 

from the figure that the along-track position estimate converges to within ±100 
feet in less than two seconds. 

The cross-track position error residual given in Figure 5.10 shows that the 
cross-track position estimate converges to within ±5 feet in less than one second. 
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Figure 5.11 illustrates the altitude error residual. It may be seen that it con- 
verges to within ±5 feet in less than six seconds. 



Figure 5.11: Altitude error using Algorithm IV. 


The along-track and cross-track velocity error residuals are shown in Figures 
5.12 and 5.13. These figures show that the along-track velocity estimate settles 
to within ±10 feet/second in less than three seconds and the cross-track velocity 
settles to within ±5 feet/second in less than a second. 

The sink rate error residual is shown in Figure 5.14. It may be observed that 
the sink rate error is reduced to within ±5 feet/second in less than three seconds. 

The position error residuals are summarized in Table 5.4. The notation and 
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units are same as those in Table 5.2. Comparing Table 5.4 to Table 3.2, it may be 


Table 5.4: Algorithm IV Results 


Category 

Lateral (j/b) 

Vertical ( h ) 

CAT I 

±1.32 

±5.77 

CAT II 

±0.69 

±2.82 

CAT Ilia 

±0.99 

±1.74 

CAT Illb & c 

±0.34 

±0.77 


noted that Algorithm IV along with the Kalman filter generates position estimates 
which meet the navigation accuracy requirements for all Categories. Moreover, 
Algorithm IV results are better than Algorithm II results with w\ — 1 and w 2 = 0 
at lower altitudes and they are far superior to the results obtained using Algorithm 
I with w-i = 1 and w 2 = 1 at all altitudes. 


5.3 Summary 

The main theme of this chapter was to exploit the structures in the model 
and in the observed lighting layouts to enable direct estimation of the position 
components. The structures formed by the threshold bar lights, left and right 
edge lights and the centerline lights were discovered by arranging the lights in a 
non-decreasing order. Two methods which utilize the correspondence between the 
structures in the model and in the observed lighting layouts were described. The 
first method required an onboard altimeter and the yaw, pitch and roll angles. 
The second method only required knowledge of the aircraft attitude angles. It 
was shown that the second algorithm is able to estimate all three position com- 
ponents without any iterative computations. A six-state Kalman filter was used 
for integrating the information derived from these methods to improve the aircraft 
position estimates and to estimate velocity components along the descent path. 




Chapter 6 


Kalman Filter Integrated 
Methods 


The algorithms described in this chapter are based on the second solution 
family described in Chapter 3. Note that the position determination algorithms in 
the two previous chapters were based on the first solution family. The reason for 
using the second solution family is that they are suitable for use as the basis for 
developing predictor /corrector methods. The techniques are based on matching 
the features between camera image and model based image. Since the model based 
image is synthesized by using the model of the airport lighting, camera model and 
the estimated aircraft position and orientation, the difference between the model 
image based and camera image based measurements can be used in a feedback 
scheme for driving the differences to zero, thereby recovering the position and 
orientation of the aircraft. This concept, coupled with the fact that the aircraft 
equations of motion temporally relate the aircraft states, allows the state estima- 
tion problem to be cast as a Kalman filtering problem. A summary of the Kalman 
filtering algorithm is given in Appendix A of this report. 

In the previous two chapters, a six-state Kalman filter was used for integrating 
the aircraft motion with the position estimates generated by the algorithms. The 
position estimates were used as inputs to the Kalman filter. In this sense, the 
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integration of the aircraft dynamics with the position determination algorithms was 
external. The Kalman filters used in this chapter attempt a much deeper level of 
integration between the vehicle dynamics and the observed image sequence. Unlike 
the algorithms in the previous chapters, correspondence between airport lighting 
layout model and the observed images is achieved by minimizing the difference 
between the camera image and the model image using the Kalman filter. 

The Kalman filters employed in the present work have their basis in the re- 
search reported in Reference [90]. In that work, three different Kalman filters were 
computed for the position determination problem. Kalman filters in the sensor 
coordinate system and in the inertial coordinate system were set up and compared 
for an image based ranging problem. It was shown that the Kalman filters in the 
sensor coordinate system and in the inertial coordinate system had comparable 
accuracies. However, the formulation in the inertial coordinate system was much 
easier to implement. An additional advantage is that since the translational states 
are linear in the inertial frame, the process update part of the Kalman filter is 
linear when formulated in this coordinate system. 

The ensuing sections describe three Kalman filtering algorithms. The first 
and second algorithms described in this chapter, Algorithm V and Algorithm VI, 
assume that the orientation angles t/>, 6 and <j> are known. These algorithms are 
designed to estimate the three runway relative position components. Algorithm V 
uses information from the camera image and the airport lighting model for position 
determination. 

Algorithm VI fuses the image and airport lighting model information with 
aircraft position estimates obtained from the Global Positioning System (GPS). 
GPS is an emerging technology for navigation. It is a satellite based navigation 
system that determines aircraft position. There is a strong movement within the 
aeronautical community to incorporate GPS receivers in every aircraft. GPS based 
position estimates can be integrated with the image based position determination 
algorithms to improve the accuracy and robustness of position estimates. Such 
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an integration can also help develop a navigation instrument that synergistically 
exploits all the information sources in an airport environment. This forms the 
motivation for the development of Algorithm VI. Algorithm VII is designed to 
provide both runway relative position and orientation. Unlike Algorithms V and 
VI, this algorithms only assumes the availability of the roll orientation angle. Ear- 
lier versions of algorithms VI and VII have been reported in [16] and [17]. 

Since several measurement equations and computations are common to all 
the algorithms, they are discussed first in this chapter. The raw data derived from 
the image consists of the image coordinates of every light visible in the image. 
Thus, the coordinates of a light p, [u p ,v p ] T , are available in the raw data. The raw 
measurements have a position uncertainty due to the errors introduced during the 
imaging process. These errors can be modeled as: 

Xl p — U p + Tju ( U * 1 ) 

V p = V p + T] v (6.2) 

where, t/ u and T] v represent pixel position uncertainty. For modeling purposes, i] u 
and r) v can be assumed to be independent scalar white noise processes. 

Clearly, the coordinates of individual lights contain little or no information 
regarding the shape or size of lighting layout. In order to incorporate this infor- 
mation, secondary measurements obtained by combining the coordinates of some 
or all the lights are needed. One of the ways of generating information regarding 
the shape and size is to assume that the image coordinates are random variables 
and construct aggregation formulae that characterize the shape and size. With 
this notion, size and shape can be related to the distribution of the two random 
variables. In that case, the characteristics of the distribution such as mean, vari- 
ance and higher-order and moments can be used to establish the size and shape 
of the object in the image. Additional measures such as the correlation coefficient 
and the eigen values of the covariance matrix or matrix singular values can also be 
used. 
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Extensive numerical experiments revealed that, the following six image-based 
measurements are found to be useful for characterizing the shape and size of the 
airport lighting layout: 

*i = £ Up/M (6.3) 

1 <p< A/ 


Z2 = E *p/M (6-4) 

\<p<M 



In these equations, M is the number of lights detected in the image; u p and v p are 
the measured coordinates of the individual light sources. The last three measure- 
ments, Equations (6.6) through (6.8) may be compactly written as: 



where, j = 4,..., 6 correspond to r = 0, t = 45 and r = 90 degrees. The 
secondary measurements in Equations (6.3) through (6.8) aggregate the size and 
shape information about the airport lighting layout. 

Physically, the first and second measurements are the arithmetic means of 
positions of the observed airport lights in the image plane. They can also be 
thought of as the coordinates of the centroid of the observed light distribution in the 
image plane. The third measurement gives the the mean distance of the observed 
lights from the origin of the image coordinate system. The fourth measurement 
is the moment about an axis parallel to the u image axis and passing through 
the centroid defined by the mean. Similarly, the fifth measurement is the moment 
about an axis passing through the cluster centroid and is inclined at 45 degrees 
to both the u and v image axes. Finally, the sixth measurement is the moment 



158 


about an axis which passes through the mean and is parallel to the v image axis. 
The fourth, fifth and sixth measurements are the square roots of the elements of 
the covariance matrix. Specifically, the fourth and sixth measurements are the 
standard deviations of the v p and u p coordinates. The fifth measurement is the 
square root of the weighted sum of the variances of u p and v p and their covariance. 

It may be noted that the six measurements z\ through z 6 in Equations (6.3) 
through (6.8) can only assume positive values. The first two measurements are 
always positive due to the choice of the origin of the image plane defined in Figure 
3.6. The other four measurements are always positive due to the use of sums of 
squares. Note that the squaring operation tends to decrease the effects of small 
numbers, while amplifying the influence of large numbers. Use of square root oper- 
ation in the secondary measurements Z 3 through z 6 prevent the position estimation 
algorithms from being biased towards large errors. 

As described in Appendix A, in order to utilize these secondary measurements 
using a Kalman filter, the measurements have to be related to the aircraft states 
as follows: 

Z = h + Cz (6.10) 

where, Z is the 6 x 1 vector of secondary measurements z\ through z 6l h is the 
6x1 state dependent measurement model vector and ( z is the measurement noise 
vector. The components of the h vector, hi through /*6, are: 


*i = £ u p /N (6.11) 

1 <p<N 

Aa « £ V P/ N (6.12) 

1<P<N 

*3 = £ (vGj + (6.13) 

1 <p<N 

K = ,/ £ (v„-/y7JV (6.14) 

y \<p<n 

h 5 = I £ (lip - hi + Vp - /i 2 ) 2 /(2iV) (6.15) 

V i<p<w 

he = .1 £ (Up-h^/N (6.16) 

y 1 < P <N 
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where N is the number of model lights within the synthesized image and u p and 
v p are the image coordinates of the lights. 

For computation of the Kalman gain and propagation of the state error covari- 
ance matrix, a linearized measurement matrix H(fc), is required. The elements of 
the H(/c) matrix can be obtained by evaluating the partial derivatives of the mea- 
surement model vector with respect to the estimated states. The partial derivatives 
used in the computation of the measurement matrix H(fc) can be obtained using 
the perspective projection Equations (3.41) and (3.42) as: 
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(6.17) 

(6.18) 

(6.19) 

( 6 . 20 ) 
( 6 . 21 ) 
( 6 . 22 ) 

(6.23) 

(6.24) 

(6.25) 

(6.26) 
(6.27) 
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dv p _ r 

a i J 


fdzcp) 

1 ( dx cv\ 

x ' c p \ d<t> ) 

1 Zc V \ d<t> ) 


(6.28) 


d<t> ~ Cp 

where the six position and orientation states are xj,, yb, Zb, Vb $ and 4>- f is the 
focal length of the camera. The partial derivatives of x cp , y Cp and z cp with respect 
to the position states are obtained from Equations (3.7) through (3.9) as: 


dx cp 

dx b 

— 1*1 

(6.29) 

II 

A 

-r 2 

(6.30) 

dx cp 

dz b 

-r 3 

(6.31) 

dy Cp 

dx b 

-r 4 

(6.32) 

dy Cp 

dyb 

~r 5 

(6.33) 

% P 

dz b 

-r 6 

(6.34) 

dz cp 

dx b 

~r 7 

(6.35) 

dz cp 

dyb 

-r$ 

(6.36) 

dz cp 

dz b 

-r$ 

(6.37) 


Similarly, the partial derivatives with respect to the orientation states are also 
obtained from Equations (3.7) through (3.9) as: 

If - <«,-*>&+<*-*>£+<*-*>£ (6 - 38 > 

d -w = +<*-*>£ +<%-*>& (« j ») 

If ~ < 6 - 40 > 

If = + + (6.41) 

If = (x ’- x>) W + iyr - yt) W + { ‘ r - Zb) lii (6 ' 42) 

~§f = + ( Vr ~ y ^~dt +< '* r ~ z ^Tw (6 ' 43) 
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d ~w = 

^-cp , sdr 7 dr 8 , , ,dr 9 

= <*■> " *»> W + ^ - »‘>a* + < 2 " - 2 ‘> H 


(6.44) 

(6.45) 

(6.46) 


Here, rj through r 9 are the components of the transformation matrix from inertial 


frame to camera frame given by Equation (3.3). It may be noted that the partial 
derivatives in Equations (6.29) through (6.46) can be computed for each light p in 
the image, synthesized from the model of the airport lighting using the estimated 


position and orientation states. 


The partial derivatives of h\, h 2 and h 8 measurements can be computed using 
Equations (6.17) through (6.28) and Equations (6.11) through (6.13) as follows: 
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(6.47) 

(6.48) 

(6.49) 

(6.50) 

(6.51) 

(6.52) 

(6.53) 

(6.54) 

(6.55) 

(6.56) 

(6.57) 



162 


d/i 2 

= - V 

<9u p 

dp 

dp 

dh 3 

- i r 

"pfe + Vvk 

dx b 

/V " 

iV i< p <n 

\/ u p + U P 

dh 3 


9u p . 

dyb 

V^p + v p 

dh 3 


8up . dv v 

U P dz h + V P dz b 

dz b 

= — 

l<p<JV 

\/ U l + V l 

dh 3 

- — E 

du p , dv p 

U V d ip + V P dip 

dip 

*1 t<N 

\j u l + v l 

dh 3 

- 1 r 

du p , dv p 

u p d0 4- v p 3e 

d0 

^1<7<JV 

\/ u l + v l 

dh 3 

- — y 

dup , dvp 

u p~df + v p ~ af 

d p 


\J u l + v l 


( 6 . 58 ) 

( 6 . 59 ) 

( 6 . 60 ) 

( 6 . 61 ) 

( 6 . 62 ) 

( 6 . 63 ) 

( 6 . 64 ) 


where, N is the number of lights. Similarly, the partial derivatives of /i 4 , h 5 and 
he can be evaluated as: 


dkj Ei< P <Ar[(w P - Msinr, + (v p - h 2 ) cos Tj][(§^ - f^)sinTj + (§jj - f^jcosTj] 


dxb < p <n[( u p - sin Tj + (v p - h^jcosr^ 

( 6 . 65 ) 

dhj Ei<p<v[(« P “ M sinT,- + (tfr - h 2 ) cos rj][(fg - g^jsinTj + (fg - f^jcosr,-] 
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( 6 . 66 ) 
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( 6 . 69 ) 
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dhj 

dcp 


Ei< P <iv[( u P - M sinr j + ( v p ~ M COST j][(?^ ~ j £) sin T i + irjf ~ jf ) cos r 
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— h\) sin Tj -f ( v p — hi) cos Tj ] 2 


(6.70) 


with Tj = 0, 45 and 90 degrees for j = 4, 5 and 6 respectively. 

The complete observation matrix H(fc) can now be written in terms of the 
partial derivatives of hi through h 6 . The details of the H(fc) matrix for each of 
the three algorithms will be described in the ensuing sections. 

In addition to the measurement model, a linearized discrete time dynamic 
model is required for implementing the Kalman filter described in Appendix A. 
All the three algorithms described in this chapter use a discrete time dynamic 
equation of the form: 

X{k + 1) = *(k)X(k) (6.71) 

X(*) is the current state vector and X(£ + 1) is the state vector at the next 
sample instant. $(k) is the state transition matrix. Since different state vectors 
are employed in each of the three algorithms, the state transition matrices are 
given separately in the following sections. 


6.1 Algorithm V 

This algorithm assumes that the aircraft yaw, pitch and roll orientations Vq 9 
and (j) are known. This algorithm estimates the three aircraft position components 
Xb, Vb, and the three aircraft velocity components Vb x , Vb y and Vb z . The 6x6 
state transition matrix is obtained by assuming that the aircraft velocity vector 
components remain constant, and that the velocity to position integration can be 
adequately approximated by the Euler integration [35] method. Thus, the state 
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transition matrix is of the form: 




1 0 0 A t 
0 10 0 
0 0 10 
0 0 0 1 
0 0 0 0 
0 0 0 0 


At is the update time step. 

The observation matrix is: 


dhx 

dh} 

a/M 

dx b 

dyb 

dz b 

dh.2 

dh 7 

dhi 

dx b 

&yb 

dz b 

dhz 

a /ia 

a /la 

dx b 

&yb 

dz b 

a/u 

a/i 4 

dhi 

dx b 

&yb 

dz b 

dh^ 

dh* 

dh>, 

dx b 

dyb 

dz b 

dh* 

a/i fi 

dhe 

Sx b 

&yb 

dz b 


0 0 
At 0 

0 At 

0 0 

1 0 
0 1 


0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 


(6.72) 


(6.73) 


In addition to the state transition matrix and the observation matrix, several 
other matrices are required for implementing the Kalman filter as described in 
Appendix A. These are defined below. 

The control input vector is a 6 x 1 null vector. For convenience, the input 
distribution matrix T(k) and the disturbance distribution matrix T^(k) are chosen 
to be 6 x 6 identity matrices. The process noise covariance matrix Q(k) is a 6 x 6 
null matrix. The dimension of the measurement noise covariance matrix R is 
6x6. In the current implementation, the diagonal elements are set to the variance 
of .25/ A t corresponding to the standard deviation of pixel position uncertainty of 
0.5 pixels. At is the measurement update time step. Note that the variance of 
pixel position uncertainty has been divided by At to convert to the discrete time 
case. 

In order to begin state estimation, the 6x1 state vector X(fc) and its error 
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covariance matrix P (A), of size 6x6, have to be initialized. The Kalman filtering 
algorithm described in Appendix A is then used for recursive state estimation. 

6.1.1 Results Using Algorithm V 

The position and velocity estimates generated by Algorithm V are described 
in this section. As before, the simulation scenario and the initial conditions were 
taken to be same as those used for the previous algorithms described in Chapters 
4 and 5. 

The along-track position error residual is shown in Figure 6.1. This figure 


o 



Figure 6.1: Along-track position error using Algorithm V. 
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shows that the along-track position estimates converge to within ±100 feet in less 
than two seconds. 

The cross-track position error residual portrayed in Figure 6.2 shows that the 
cross-track position estimate converges to within ±5 feet in less than one second. 



Figure 6.2: Cross-track position error using Algorithm V. 

Figure 6.3 illustrates that the altitude estimates converge to within ±5 feet 
in less than two seconds. 

The along-track velocity error residual is shown in Figure 6.4. It may be seen 
that the along-track velocity estimates settle to within ±10 feet/second in less 
than six seconds. 

The cross-track velocity error residual in Figure 6.5 shows that the cross-track 
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velocity settles to within ±5 feet/second in less than a second. 



Figure 6.5: Cross-track velocity error using Algorithm V. 

The sink rate error residual in Figure 6.6 shows that the sink rate error also 
is reduced to within ±5 feet/second in less than a second. 

The position error residuals are summarized in Table 6.1. The notation and 
units are same as those in the previous tables. Comparing Table 6.1 to Table 3.2, 
it may be seen that the position estimates resulting from Algorithm V meet the 
navigation accuracy requirements for all Categories. Table 6.1 shows that amongst 
the five algorithms discussed so far in this report, Algorithm V provides the most 
accurate aircraft position estimates. 



Time ( s ) 

Figure 6.6: Sink rate error using Algorithm V. 


Table 6.1: Algorithm V Results 
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6.2 Algorithm VI 

Like the previous algorithm, this algorithm also assumes that the aircraft 
orientation angles are known. In addition to the six measurements described in 
Equations (6.3) through (6.8), this algorithm assumes that three position compo- 
nents are provided by an onboard Global Positioning System (GPS) receiver. 

The reasons for integrating the GPS with the vision based position determi- 
nation algorithm are as follows. The GPS provided position of the aircraft can 
be used to initialize the Kalman filter. If during descent, airport lights are cut off 
due to foreground occlusion, the integrated algorithm would continue to provide 
estimates of the aircraft position using GPS measurements. If the GPS signals 
are blocked due to terrain obstacles such as mountains or buildings, the integrated 
system would continue to estimate aircraft position using the vision based system. 
Since commercial GPS systems with standard position service have an accuracy 
of 325 feet horizontally 95 percent of the time [30] and 560 feet vertically, the 
bias and the noise in the GPS position can be reduced by the integrated system. 
Thus, integration of the vision system with the GPS is motivated by robustness. 
Moreover, the integrated navigation system synergistically exploits the available 
data sources for position estimation. 

Since the location of the airport is known, the GPS-based aircraft position 
can be used for estimating the runway relative aircraft position components x*,, y 
and Zb. These can be modeled as: 


%b — ^6 1 b x "t" Vx 

(6.74) 

Vb = Vb + by + Vy 

(6.75) 

z b = z b + b z + tj z 

(6.76) 

The GPS measurement model includes the bias terms b x , 

by and 6 Z , and white 


noise terms, r \ x , r\ y and rj z . The GPS-based position components are directly used 
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as three additional measurements z 7 , z 8 and z$ in the present formulation. Thus, 


z 7 

= x b 

(6.77) 

Z S 

= Vb 

(6.78) 

z 9 

= h 

(6.79) 


Based on Equations (6.74) through (6.76), the components of the measurement 
model vector /i 7 , hg and are: 


h 7 

— Xi + b x 

(6.80) 

hs 

= 2/6 + by 

(6.81) 

hg 

- z h + b z 

(6.82) 


Assuming the bias terms to be constant during approach and landing, the 
following state equations can be used to model the bias terms: 


b x {k + 1) = b x (k) (6.83) 

b y (k + 1) = b y (k) (6.84) 

b z {k + 1) = b z (k) (6.85) 

The state vector for the algorithm consists of aircraft position and velocity 
components, together with the bias vector in the GPS. Thus, 

X = [Xb,Vb,B] r ( 6 . 86 ) 


is the aircraft position vector, is the aircraft inertial velocity vector and B 
is the GPS bias vector. Assuming constant inertial velocity and Euler integration 
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for state vector propagation results in the following 9x9 state transition matrix: 


*(*) = 


1 0 0 At 


0 1 0 
0 0 1 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 


0 

0 

1 

0 

0 

0 

0 

0 


0 

At 

0 

0 

1 

0 

0 

0 

0 


0 0 0 0 
0 0 0 0 
At 0 0 0 
0 0 0 0 
0 0 0 0 
10 0 0 
0 10 0 
0 0 10 
0 0 0 1 


(6.87) 


At is the update time step. 

With 9 elements of the measurement model vector, six image-based mea- 
surements defined in Equations (6.11) through (6.16) and three GPS positions in 
Equations (6.80) through (6.82), the complete observation matrix is: 


H(*) = 


dhy 

dx b 

dhi 

dVb 

dhx 

dz b 
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0 

0 

dh 7 

dx b 
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0 

0 
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0 
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dyb 
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0 

0 

0 

0 

0 

a/i fi 

dx b 

dhfi 

dyb 

dhe 

dz b 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 


( 6 . 88 ) 


The partial derivatives of hi through with respect to position components were 
discussed earlier in Equations (6.47) through (6.70). 

The implementation of this Kalman filter differs significantly from Algorithm 
V because the image-based and GPS-based measurements are available at two 
different rates. The image-based measurements are typically available at the rate 
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of ten times a second, while the GPS measurements are available once a second. 
Hence, during one second, six image-based measurements Z\ through z$ are avail- 
able ten times, and the three GPS-based measurements z-j through zg are available 
once. A multi-rate formulation of the Kalman filter [63] is required to deal with 
the changing dimension of the measurement vector. 

The multi-rate formulation requires the observation matrix H(£) and the 
measurement noise covariance matrix R(&) to change with the number of mea- 
surements available at any one measurement epoch. For instance, the last three 
rows of the H(fc) matrix in Equation (6.88) are eliminated when GPS measure- 
ments are unavailable. The noise covariance matrix R(fc) is a 9 x 9 diagonal matrix 
when both GPS-based and image-based measurements are available. The first six 
diagonal elements contain the pixel noise variances, while the remaining three diag- 
onal elements contain GPS measurement noise variances. The standard deviation 
of the pixel noise is assumed to be 0.5 pixels and the GPS measurement noise is 
three meters. It may be noted that the pixel variance is scaled by 0.1 second while 
the GPS measurement variance is scaled by 1 second. When only the image-based 
measurements are available, the last three rows and columns of the R(fc) are elim- 
inated in order to reduce its dimension to 6. The remaining matrices required for 
implementing the Kalman filter are defined below. 

The control input vector is a 9 x 1 null vector. The input distribution matrix 
r(fc) is set to be a 9 x 9 identity matrix. The disturbance distribution matrix 
r<i(fc) is a 9 x 9 identity matrix. The process noise covariance matrix Q(A:) is a 
9x9 diagonal matrix. The diagonal elements are chosen as “tuning parameters” 
for the Kalman filter [35]. 

In order to begin the state estimation process, the 9x1 state vector X(£) and 
its error covariance matrix P (k), of size 9x9, have to be initialized. Following 
the standard procedure in Kalman filtering, the state error covariance matrix is 
initialized by placing large variance values in the diagonal locations and setting 
the off-diagonal terms to zero. The Kalman filtering algorithm in Appendix A is 
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then used for recursive state estimation. 

On a final note, it may be verified from Equation (A. 5) in Appendix A that 
the dimension of the Kalman gain matrix changes from 9 x 6 to 9 x 9 depending 
on whether image-based measurements and/or GPS measurements are available. 
Thus, the structure of the Kalman filter is changed to accommodate the change 
in the number of measurements. This way, the Kalman filter is always updated at 
the fastest measurement rate. 

6.2.1 Results using Algorithm VI 

The position, velocity and GPS bias estimates obtained by Algorithm VI are 
described in this section. The estimates were obtained for the same simulation 
scenario used in all the previous algorithms. The GPS position estimates were 
used for initializing the Kalman filter. In order to simulate GPS measurements, 
position bias with a uncertainty of ±100 meters and Gaussian white noise with 
a standard deviation of three meters were added to the aircraft position vector. 
The GPS position bias errors of 82 feet in the along-track position, 177 feet in the 
cross-track position and 142 feet in the altitude were assumed for generating the 
results given in this section. 

The along-track position error residual is shown in Figure 6.7. This figure 
shows that the along-track position estimates converge to within ±100 feet in two 
seconds. 

The cross-track position error residual given in Figure 6.8 shows that the 
cross-track position estimate converges to within ±5 feet in less than one second. 

Figure 6.9 illustrates that the altitude estimates converge to within ±5 feet 
in less than two seconds. 

The along-track velocity error residual is shown in Figure 6.10. It may be 
seen that the along-track velocity estimates settle to within ±10 feet/second in 
less than two seconds. 

The cross-track velocity error residual in Figure 6.11 shows that the cross- 
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track velocity settles to within ±5 feet/second in less than a second. 



Figure 6.11: Cross-track velocity error using Algorithm VI. 

The sink rate error residual shown in Figure 6.12 reveals that the sink rate 
error also is reduced to within ±5 feet/second in less than two seconds. 

The GPS bias error residuals in the along-track and cross-track positions and 
in the altitude are shown in Figures (6.13), (6.14) and (6.15). The GPS bias 
estimates in the along-track position settles to within ±20 feet in less than six 
seconds. The GPS bias estimates in the cross-track position and altitude converge 
to within ±10 feet in less than 17 seconds and three seconds respectively. The 
altitude bias error increases beyond ±10 feet when the aircraft is 36 feet above the 


runway. 
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The position error residuals are summarized in Table 6.2. The notation and 
units are same as those in the previous tables. The usefulness of the integrated 


Table 6.2: Algorithm VI Results 


Category 

Lateral ( y &) 

Vertical ( h ) 

CAT I 

±0.02 


CAT II 

±0.01 


CAT Ilia 



CAT Illb & c 

±0.02 



navigation system in aircraft operations can be assessed by comparing the achieved 
accuracy with the desired accuracies listed in Table 3.2. It may be observed from 
Table 6.2 that the present navigation scheme meets the objectives for all the Cat- 
egories listed in Table 3.2. Table 6.2 shows that Algorithm VI is able to provides 
accurate aircraft position estimates. Moreover, Algorithm VI improves the GPS 
bias estimates considerably. 

The results presented here show that the integrated GPS and image based 
algorithm is as accurate as the pure image based algorithm. The accuracy of the 
integrated algorithm can be further improved with a GPS receiver with a faster 
update rate. The integrated algorithm is expected to be fault tolerant in the event 
of imaging system failure. 


6.3 Algorithm VII 

In the two previous algorithms, the orientation angles were assumed to be 
known. In this final algorithm of the present work, the orientation angles are 
included as additional states to be estimated. The objective is to examine the 
degree to which the image based position determination concept can be extended. 
As discussed in Section 2.3. 7.1, it is unlikely that all three attitudes can be reliably 
estimated without explicit correspondence between several points in the camera- 
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based and model-based images. Hence, in the present work, the attempt will be 
to estimate the aircraft pitch and yaw angles, with the assumption that the roll 
angle is available from measurements. It needs to be mentioned here that the 
attempts at including all the three attitudes did not yield useful results. With the 
addition of orientation states, this algorithm can be considered to be an extension 
of Algorithm V. 

Using Euler Equations (3.24) through (3.26) describing the attitude kinemat- 
ics, the discrete time state equation for attitudes can be obtained as: 


( rP(k+ 1) ^ 


1 0 cos <p(k) sec0(k)At 

sin <^(A:) sec 6(k)At 


' ip(k) ^ 

0{k + 1) 


0 1 — sin<^(A;)Af 

cos (p(k)At 


0(k) 

r(k + 1) 


0 0 1 

0 


r(k) 

^ <l{k +1) J 


r 

O 

O 

o 

1 


v / 


(6.89) 

where r and q are the yaw and pitch angular rates, and At is the time step. The yaw 
and pitch angular rates are assumed to be constant during the landing phase. The 
Euler equation shown in Equation (6.89) was discretized using forward differences 
which resulted in an explicit scheme. Note that the Euler equation can also be 
discretized using backward differences to yield an implicit scheme. The advantage 
of an implicit scheme is that it is unconditionally stable with respect to step size 
[45]. An explicit scheme becomes unstable for large step sizes. The disadvantage 
of an implicit scheme is that it requires several computations at each time step. 
Due to this reason, the explicit scheme has been used here. 

Since the state vector is formed by combining the translational position and 
velocity states xj, y*,, Zb, Vb x , Vb y and Vb z with the rotational states: ip, 0, r and q , 
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the 10 x 10 state transition matrix is: 

100 At 0000 0 0 

0100 At 000 0 0 

00100 At 00 0 0 

00010000 0 0 

00001000 0 0 

00000100 0 0 

0000 0 010 cos 4>(k) secO(k)At sin<i>(k)secd(k)At 

0000 0 001 —sin<f>(k)At cos 4>(k) At 

00000000 1 0 

00000000 0 1 

(6.90) 

Since the quasi-steady state approximation in the orientation angles influence the 
state transition matrix, a further approximation is introduced in the propagation 
of the error covariance matrix during the time update step described in Equation 
(A. 8) in Appendix A. 

Using the measurement Equations (6.11) through (6.16) and the ten states, 
the 6 x 10 observation matrix can be constructed as: 
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(6.91) 


Other vectors and matrices needed 


for implementation of the Kalman filter 


are as follows. The control input vector is a 10 x 1 null vector. The input dis- 
tribution matrix T(k) is a 10 x 10 identity matrix. The disturbance distribution 
matrix Td(A:) is also a 10 x 10 identity matrix. The dimension of the process noise 
covariance matrix Q (k) is 10 x 10. 
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Just as in algorithm V, six measurements z\ through ze defined in Equations 
(6.3) through (6.8) are used. Thus, the 6x6 measurement noise covariance matrix 
R (k) in Algorithm V is also employed here. 

The lOx 1 state vector X(Ar) and its error covariance matrix P(Ar), of size 10x10 
are initialized using an approach identical to that in the six previous algorithms. 
The Kalman filtering algorithm is then used for recursive state estimation. 

6.3.1 Results using Algorithm VII 

The estimates of position, velocity, yaw and pitch orientation angles, and yaw 
and pitch body rates generated by Algorithm VII are described in this section. The 
estimates are obtained for the same simulation scenario used in all the previous 
algorithms except that the yaw orientation angle is set to -10 degrees. As before, 
the pitch orientation angle is set to -3 degrees. Errors of 1000 feet in the along- 
track position Xf, and 100 feet in the cross-track position yj, and in the altitude 
— Zb are assumed for initializing the Kalman filter. Initial values of the velocity 
components, the yaw and pitch orientation angles, and the yaw and pitch body 
rates are set to zero. 

The along-track position error residual is shown in Figure 6.16. It may be 
observed that the along-track position estimates converge to within ±100 feet in 
five seconds. 

The cross-track position error residual is shown in Figure 6.17 reveals that the 
cross-track position estimate converges to within ±5 feet in less than one second. 

Figure 6.18 illustrates that the altitude estimates converge to within ±12 feet 
in about six seconds. 

The along-track velocity error residual is shown in Figure 6.19. It may be 
observed that the along-track velocity estimates settle to within ±10 feet/second 
in 11 seconds. 

The cross-track velocity error residual in Figure 6.20 shows that the cross- 
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track velocity settles to within ±5 feet/second in about a second. 



Figure 6.20: Cross-track velocity error using Algorithm VII. 

The sink rate error is reduced to within ±5 feet/second in about five seconds 
as can be observed in Figure 6.21. 

The aircraft yaw and pitch attitude error residuals are presented in Figure 
6.22 and 6.23 respectively. Both these errors settle to within ±0.2 degree in less 
than a second. 

Figures 6.24 and 6.25 show that the yaw and pitch body rate error residuals 
settle to within ±0.1 degrees/second in under two seconds. 

The position error residuals are summarized in Table 6.3. The notation and 
units are same as those in the previous tables. By comparing Table 6.3 to Table 
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Figure 6.21: Sink rate error using Algorithm VII. 


Table 6.3: Algorithm VII Results 


Category 

Lateral (j/j) 

Vertical ( h ) 

CAT I 

±0.60 

±11.29 

CAT II 

±0.25 

±10.10 

CAT Ilia 

±0.85 

±2.54 

CAT Illb & c 

±1.36 

±0.86 
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3.2, it may be seen that the results obtained using Algorithm VII only satisfy 
Category I landing requirements. The results presented in this section show that 
Algorithm VII provides very accurate estimates of the aircraft attitude angles. 


6.4 Summary 

Three Kalman filtering centered algorithms were presented in this chapter. 
The first algorithm was designed to provide estimates of the aircraft position and 
velocity components. The second algorithm was designed to provide estimates of 
the aircraft position, velocity and GPS position bias components. Finally, the third 
algorithm was designed to provide estimates of the aircraft position and velocity 
components, yaw and pitch orientation angles, and yaw and pitch body rates. 

All three algorithms used six shape features of the airport lighting layout 
as image-based measurements for the Kalman filters. Additionally, Algorithm 

VI used the three aircraft position components provided by the GPS receiver. 
Except for Algorithm VII, the other two algorithms also needed yaw, pitch and 
roll orientation angles. Algorithm VII only required the roll orientation angle. 

Results were obtained to demonstrate the performance of the Kalman fil- 
ters. It was shown that Algorithms V and VI are able to provide aircraft position 
estimates which meet Category I, II, Ilia, Illb and IIIc navigation accuracy re- 
quirements. Aircraft position estimates generated using Algorithm VII were only 
able to meet Category I navigation accuracy requirements. However, Algorithm 

VII was able to provide highly accurate estimates of the yaw and pitch orientation 
angles. 



Chapter 7 


Contributions of the Report and 
Future Work 


The contributions of this report and future research directions are discussed 
in this chapter. 

7.1 Contributions of This Report 

This report has explored the development of machine vision based pilot aids 
to help reduce night approach and landing accidents. The research focus was on 
developing an onboard instrument that complements the existing cockpit instru- 
mentation. 

The techniques developed during the course of this research were motivated 
by the desire to use the existing information sources to derive more precise aircraft 
position and orientation information. During night landing, the information source 
used by pilots for obtaining aircraft position and orientation information is the 
airport lighting layout. The fact that airport lighting geometry is known and since 
the images of the airport lighting can be acquired from an onboard camera, machine 
vision technology can be used for synthesizing a landing aid. Use of a machine 
vision system has several advantages. Firstly, such systems are not susceptible to 
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optical illusions. Moreover since the camera is a passive imaging device, it does not 
cause interference with the ground based equipment or with equipment onboard 
other aircraft. Finally, lowering costs of electro-optical cameras and real-time 
computer systems have made this technology attractive. Even if Global Positioning 
System (GPS) receivers become cheaper and more accurate, an integrated machine 
vision and GPS system would be a much more robust landing aid. The machine 
vision based system could also serve as a back-up landing system. 

The main contribution of this research are the synthesis of seven navigation 
algorithms based on two broad families of solutions. The first family of solution 
methods comprise of techniques that reconstruct the airport lighting layout from 
the camera image and then estimate the aircraft position components by compar- 
ing the reconstructed lighting layout with the known model of the airport lighting 
layout. The second family of methods consist of techniques that synthesize the 
image of the airport lighting using a camera model and the known model of the 
airport lighting layout and then estimate the aircraft position components by com- 
paring this synthesized image with the actual image of the airport lighting acquired 
by the onboard camera. 

Algorithms I through IV belong to the first family of solutions while Algo- 
rithms V through VII belong to the second family of solutions. These algorithms 
can further be classified as parameter optimization methods, feature correspon- 
dence methods and Kalman filter centered methods respectively. Algorithms I 
and II are parameter optimization methods. Algorithms III and IV are feature 
correspondence methods. Algorithms V, VI and VII are Kalman filter centered 
methods. Figure 7.1 summarizes the algorithm classification. 

Figure 7.1 shows the two classes of machine vision based landing aid developed 
in this report. First class provide only position information and second category 
provide both position and orientation information. Algorithms I. and III provide 
the aircraft x and y inertial position components. They assume that the altitude 
information is available from an onboard altimeter. Algorithms II, IV, V and VI 



202 



Figure 7.1: Classification of the seven algorithms developed in this report. 
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compute all three aircraft inertial position components. Additionally, Algorithms 
V and VI provide estimates of the aircraft inertial velocity. Since Algorithm VI 
integrates image-based measurements with the position measurements from a GPS 
receiver, it also provides estimates of GPS position bias components. These can 
be used to improve the accuracy of GPS measurements. 

Algorithms I through VI all assume that the yaw, pitch and roll attitude angles 
are available. Algorithm VII provides estimates of all three runway relative aircraft 
position and velocity components, the yaw and pitch orientation angles, and the 
yaw and pitch body rates. This algorithm only assumes that the roll attitude 
angle is available. Table 7.1 summarizes the aircraft states estimated by the seven 
algorithms developed in this report. The estimated quantities are marked with 
bullet. Note that velocity estimates for Algorithms I, II, III and IV are obtained 
by using a six-state Kalman filter driven by these algorithms. In Table 7.1 x b , 
y b and z b are the inertial position components, v bx , v by and v bz are the inertial 
velocity components, 6 r , b y and b. are the GPS position bias components, ip and 0 
are the yaw and pitch orientation angles, and r and q are the yaw and pitch body 
rates. Table 7.1 also indicates the nature of the computations required for each of 
these algorithms. Algorithms I, III and IV are direct computational schemes that 
do not require iterative computations. Algorithm II, V, VI and VII use iterative 
computational schemes. 

In order to take advantage of the aircraft dynamics and the multiple images 
available along the glide path, the estimates provided by Algorithms I, II, III 
and IV were used for driving a six-state Kalman filter for providing estimates of 
the aircraft position and inertial velocity components. Algorithms V, VI and VII 
are Kalman filter centered algorithms and were designed to implicitly utilize the 
aircraft dynamics and multiple images available along the glide path. 

Results were presented to demonstrate the performance of all the seven al- 
gorithms developed in this report. It was shown that all the algorithms are able 
to meet some or all of the Federal Aviation Administration specified navigation 
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Table 7.1: States Estimated by the Seven Algorithms 
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Alg. I 

Alg. 11 
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accuracy requirements for various landing categories. Table 7.2 summarizes the 
performance of the seven algorithms in meeting the navigation accuracy require- 
ments for various FAA categories (CAT) of landing in Table 3.2. 


Table 7.2: Performance Of The Seven Algorithms 


CAT 

Alg. 1 

Alg. II 

Alg. Ill 

Alg. IV 

Alg. V 

Alg. VI 

Alg. VII 

I 

• 

• 

• 

• 

• 

• 

• 

II 

• 

• 

• 

• 

• 

• 


Ilia 

• 


• 

• 

• 

• 



• 


• 

• 

• 

• 



7.2 Practical Considerations 

The algorithms reported in this research have been validated using simulated 
image sequences. By comparing an actual image in Figure 2.2 with the simulated 
image in Figure 3.10 it may be observed that the only difference between the two 
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images is that lights in the actual image occupy regions in the image while lights 
in the simulated image appear as point sources. Thus, the actual image can be 
processed through low-level algorithms to transform it to appear similar to the 
simulated image. One may conclude that for Algorithms I, II, III and VI recon- 
struction of the airport lighting layout based on an actual image would result in 
lights occupying regions on the plane of the runway. To convert these regions to 
points, the centroids of the regions have to be found. A circular template of the 
physical dimensions of an airport light can be used to determine the centroids. 
Once the centroids are found for an initial image, the position and velocity esti- 
mates provided by the Kalman filters can be used to aid local template matching 
in subsequent images. For Algorithms V, VI and VII, the image constructed from 
the model of the airport lighting using propagated position estimates can be used 
for aiding the search for centroids of the lights detected in the actual image. Tem- 
plate matching schemes or local clustering schemes can be used for determining 
the centroids of lights in the actual image. 

The computational requirements for all the algorithms are modest. Of the 
several Kalman filters developed in this report, the largest one is a ten-state filter 
with six measurements used in Algorithm VII. Experience with ranging algorithms 
which track three position components of several hundred objects in the image 
using three-state Kalman filters for every object has shown that these algorithms 
can be made to work in real-time using inexpensive hardware [84], 

It may be noted that Algorithms II, V, VI and VII need initial values of the 
aircraft position components. In the case of Algorithm VI, the initial estimates 
are provided by an onboard GPS receiver. Algorithm I or III can be used with a 
barometric altimeter for initializing Algorithm II and Algorithm IV can be used 
for initializing Algorithms V through VII. 
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7.3 Future Work 

Since the seven algorithms developed in this report have been verified only in 
simulation, the next logical step would be to verify the performance of these algo- 
rithms using actual images of the airport lighting layout obtained by an onboard 
camera. 

This report has employed four lighting structures, the left and right runway 
edge lights, centerline lights and the threshold bar lights, in Algorithm IV. Tech- 
niques need to be developed for detecting other lighting structures in order to 
extend Algorithm IV to estimate the yaw, pitch and roll orientation angles. 

It may be possible to improve the estimation accuracy of Algorithms V, VI and 
VII by extending them to iterated Kalman filtering algorithms. This should spe- 
cially be investigated for Algorithm VII in order to improve its altitude estimation 
accuracy. 

The six features used in Algorithm VII were found to have very low sensitiv- 
ity to the roll orientation angle. Other features based on higher order moments 
should be investigated for possible estimation of the roll orientation angle. Shape 
features such as perimeter, area, eccentricity and thinness [27, 53] should also be 
investigated for improving the robustness of Algorithms V, VI and VII and for 
possible roll orientation angle estimation using Algorithm VII. 

The focus of this research was to develop a pilot aid for flight on or near the 
glide slope. Hence, it was assumed the aircraft is headed in the approach direction 
of the runway and that only the airport lights are visible in the image. This report 
has not addressed the question of initial acquisition of the airport lighting layout 
when the aircraft is not lined up in the approach direction of the runway. Since 
the heading of the airport is known from published charts [54] and the heading of 
the airplane is known from cockpit instruments such as a gyro compass, the pilots 
may be able to identify the airport. Once they identify the airport and line up 
the aircraft in the approach direction present algorithms can be initiated. It may 
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be possible to develop a pilot aid for detecting the airport by using the heading 
information and images of the ground lighting layout. Color images may also be 
beneficial for airport detection because both approach and runway lights are color 
coded. 
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Appendix A 


Kalman Filtering Algorithm 


The Kalman Filter development is based on a linear discrete time dynamical 
model of the form: 

X(k + 1 ) = *(k)X(k) + r(*)u(*) + r d (*Kx(*) (A.i) 

X(k) is the state vector, U(/c) is the control input vector, Cx(^) is a vector of 
discrete time white noise sequences with covariance Q (k) representing the process 
noise, <&(k) is the state transition matrix [35], r(fc) is the input distribution ma- 
trix and ru(&) is the process noise distribution matrix. The Measurement vector 
equation is given by: 

Z(k) = h(X(fc)) + Cz(fc) (A. 2) 

Here, Z (k) is the measurement vector, h (X(k)) is the vector of nonlinear measure- 
ment functions and £ z (fc) is measurement noise vector with covariance R(fc). Note 
that ( z (k) is assumed to be a vector of white noise sequences. 

The Kalman Filter [3] is a computational algorithm for computing optimal 
state estimates X(k) using the linear discrete time dynamical model and the mea- 
surement equations. The Kalman filter is optimal in the sense of generating un- 
biased minimum variance estimates. The filter continuously generates the state 
estimate error covariance matrix P (k). The Kalman Filter consists of two steps: 
measurement update, which improves the state estimate based on the new mea- 
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surements, and process update, which propagates the state estimate according to 
the dynamical equations. Before every measurement update step, an estimate of 
the state X(k), state error covariance matrix P (£), process noise covariance ma- 
trix Q (k) and measurement noise covariance matrix R (k) are known. The new 
measurements are used for improving the state estimate and its error covariance 
as: 


X{k) = X{k) + K(k)[Z(k)-h(X(k))] (A. 3) 

P (k) = [I — K(k)H(k)]P(k) (A. 4) 

where I is the identity matrix, K(&) is the Kalman gain matrix computed using 

K(k) = P(k)H T {k){H(k)P(k)H T (k) + R(fc)]- 1 (A. 5) 

H(fc) is the matrix of partial derivatives, representing the linear approximation to 
the nonlinear measurement functions, computed as: 

H(fc) = 3h(X)/dX| x=i (A. 6) 

The process update part of the Kalman Filter accounts for system dynamics and 
propagates the state and its error covariance until the next measurement is ob- 
tained. The propagated values are: 

X(k + 1) = *(k)X(k) + r(*)U(Jfc) (A. 7) 

P(* + l) = *(k)P(k)* T (k) + r d (k)Q(k)P d T (k) (A. 8) 

The steps A. 3 through A. 8 form the core of the Kalman filter. This algorithm is 
summarized in Table A.l. 

The sequence of steps given in Table A.l assume that the measurement and 
process updates are carried out at the same rate. The extension to the case of 
measurement update time step being an integer multiple of the process update 
time step is straight forward [35]. The procedure is more complicated when the 



222 


Table A.l: Summary of Kalman Filter Algorithm 


1. Set k = 1. 

2. Initialize X(fc), P(fc), Q(A:) and R(&). 

3. Compute h(X(k)). 

4. Compute H(fc) using Equation (A. 6). 

5. Compute Kalman gain K(fc) using Equation (A. 5). 

6. Compute X(fc) using Equation (A. 3). 

7. Compute P (k) using Equation (A. 4). 

8. Compute X(& + 1) using Equation (A. 7). 

9. Compute P(k-f 1) using Equation (A. 8). 

10. Increment k = k + 1. 


11. Return to step 3. 
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measurement update is done asynchronously. In this case the measurement arrives 
within the process update time step. One of the ways of dealing with this situation 
is to split the process update into two steps: one from the previous process update 
time to the measurement epoch and the other from the measurement epoch to the 
next scheduled process update time. This ensures that the state estimate and its 
error covariance are available synchronously. Further details of multi-rate Kalman 
filter implementation can be found in Reference [3]. 



Appendix B 


Matrices Using Aircraft 
Kinematic Models 


The matrices required for estimating the position and velocity of the aircraft 
using Algorithms I, II, III and IV along with the Kalman filter in Appendix A are 
described in this Appendix. Algorithms I and II are described in Chapter 4 and 
Algorithms III and IV are described in Chapter 5. The outputs generated by these 
algorithms are used as measurements for the Kalman filter described in Appendix 
A. 


The position and velocity estimation problem can be stated as follows. Given 
noisy measurements of the aircraft position, estimate its position components Xb, 
yb, and Zb, and its inertial velocity components Vb x , Vb y , and Vb z . 

The discrete time state transition matrix with position and velocity compo- 
nents as states can be found to be: 


*(*) = 


1 0 0 At 0 0 

0 1 0 0 At 0 

0 0 1 0 0 At 

0 0 0 1 0 0 

0 0 0 0 1 0 

0 0 0 0 0 1 


(B.l) 
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This development assumes that the velocity states are integrated using the Euler 
integration formula. Here, At is the update time step. The control input vector 
is a 6 x 1 null vector. For convenience, the input distribution matrix T(k) can be 
chosen to be a 6 x 6 identity matrix. The process noise distribution matrix 
is a 6 x 6 identity matrix. The process noise covariance matrix Q(A?) is a 6 x 6 
diagonal matrix. The diagonal elements of Q(&) are chosen as “tuning parameters” 
for the Kalman filter [35], 

The 3x1 measurement vector Z (k) consists of the components of the aircraft 
inertial position vector, y & and z\>, determined using Algorithms I, II, III and 
IV. With these measurements, position and velocity components as states, the 
measurement model matrix is: 


H(*) = 


1 0 0 
0 1 0 
0 0 1 


0 0 0 
0 0 0 
0 0 0 


(B-2) 


Since three measurements are used, the dimension of the measurement noise co- 
variance matrix R (k) is 3 x 3. 


The 6x1 state vector X(&) and its 6 x 6 error covariance matrix P(k) have 
to be initialized to begin state estimation process. The Kalman filtering algorithm 
described in Appendix A can then be used for recursive state estimation. 



